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Abstract. We complete previous investigations on the thermodynamics of self-gravitating systems by studying the 
grand canonical, grand microcanonical and isobaric ensembles. We also discuss the stability of polytropic spheres 
in connexion with a generalized thermodynamical approach proposed by Tsallis. We determine in each case the 
onset of gravitational instability by analytical methods and graphical constructions in the Milne plane. We also 
discuss the relation between dynamical and thermodynamical stability of stellar systems and gaseous spheres. 
Our study provides an aesthetic and simple approach to this otherwise complicated subject. 
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1. Introduction 

The statistical mechanics of systems interacting via long- 
range forces exhibits peculiar features such as negative 
specific heats, inequivalence of statistical ensembles and 
phase transitions. These curious behaviours have been 
first discussed in the astrophysical literature during the 
elaboration of a thermodynamics for stars (Eddington 
1926), globular clusters (Lynden-Bell & Wood 1968), 
black holes (Hawking 1974) and galaxies (Padmanabhan 
1990). They have been recently rediscovered in differ- 
ent fields of physics such as nuclear physics, plasma 
physics, Bose-Einstein condensates, atomic clusters, two- 
dimensional turbulence... The main challenge is repre- 
sented by the construction of a thermodynamic treatment 
of systems with long-range forces and by the analogies 
and differences among the numerous domains of applica- 
tion (see Chavanis 2002e and other contributions in that 
book). 

Gravity provides a fundamental example of unshielded 
long-range interaction for which ideas of statistical me- 
chanics and thermodynamics can be developed and tested. 
For systems with long-range interactions, the mean-field 
approximation is known to be exact in a suitable ther- 
modynamic limit. Therefore, the structure and stabil- 
ity of self-gravitating systems at statistical equilibrium 
can be analyzed in terms of the maximization of a ther- 
modynamical potential. This thermodynamical approach 
leads to isothermal configurations which have been stud- 
ied for a long time in the context of stellar structure 
(Chandrasekhar 1942) and galactic structure (Binney & 
Tremaine 1987). As is well-known, isothermal spheres have 
infinite mass so that the system must be confined within 



a box (Antonov problem) in order to prevent evaporation 
and make the thermodynamical approach rigorous. It is 
also well-known that isothermal configurations only corre- 
spond to metastable equilibrium states (i.e., local maxima 
of the thermodynamical potential), not true equilibrium 
states. These metastable equilibrium states are expected 
to be relevant, however, for the timescales contemplated 
in astrophysics. In particular, globular clusters described 
by Michie-King models are probably in such metastable 
states. 

The series of equilibria of finite isothermal spheres can 
be parametrized by the density contrast between the cen- 
ter and the boundary of the system. For sufficiently low 
density contrasts, the system is thermodynamically sta- 
ble. However, instability occurs at sufficiently large con- 
centrations: at some point in the series of equilibria, the 
solutions cease to be local maxima of the thermodynami- 
cal potential and become unstable saddle points. The cru- 
cial point to realize is that the onset of instability de- 
pends on the statistical ensemble considered: microcanon- 
ical (MCE), canonical (CE) or grand canonical (GCE). 
This contrasts with ordinary systems, with short range 
interactions, for which the statistical ensembles are equiv- 
alent at the thermodynamic limit, except near a phase 
transition. This suggests that gravitating systems in virial 
equilibrium are similar to normal (extensive) systems at 
the verge of a phase transition (Padmanabhan 1990). 

The thermodynamical stability of self-gravitating sys- 
tems can be studied by different technics. Horwitz & Katz 
(1978) use a field theory and write the density of states, 
the partition function and the grand partition function in 
MCE, CE and GCE respectively as a path integral for a 
formal field </). In the mean-field approximation, the inte- 
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gral is dominated by the distribution 0o which maximizes 
a specific action In the grand canonical ensemble, 

^GC£;[0] is the Liouville action. Horwitz & Katz (1978) 
solve the problem numerically and find that the series of 
equihbria become unstable for a density contrast 1.58, 32.1 
and 709 in GCE, CE and MCE respectively. It is normal 
that the critical density contrast (hence the stability of 
the system) increases when more and more constraints are 
added on the system (conservation of mass in CE, conser- 
vation of mass and energy in MCE). This field theory has 
been rediscussed recently by de Vega & Sanchez (2002) 
who confirmed previous results and proposed interesting 
developements. 

The thermodynamical stability of self-gravitating sys- 
tems can also be settled by studying whether an isother- 
mal sphere is a maximum or a saddle point of an appropri- 
ate thermodynamical potential: the entropy in MCE, the 
free energy in CE and the grand potential in GCE. The 
change of stability can be determined very easily from the 
topology of the equilibrium phase diagram by using the 
turning point criterion of Katz (1978) who has extended 
Poincare's theory on linear series of equilibria (see also 
Lynden-Bell & Wood 1968). This method is very powerful 
but it does not provide the form of the perturbation pro- 
file that triggers the instability. This perturbation profile 
can be obtained by computing explicitly the second order 
variations of the thermodynamical potential and reducing 
the problem of stability to the study of an eigenvalue equa- 
tion. This study was first performed by Antonov (1962) in 
MCE and revisited by Padmanabhan (1989) with a sim- 
pler mathematical treatment. This analysis was extended 
in CE by Chavanis (2002a) who showed in addition the 
equivalence between thermodynamical stability and dy- 
namical stability with respect to Navier-Stokes equations 
(Jeans problem). Remarkably, this stability analysis can 
be performed analytically by using simple graphical con- 
structions in the Milne plane. 

In the present paper, we propose to extend these an- 
alytical methods to more general situations in order to 
provide a complete description of the thermodynamics of 
spherical self-gravitating systems. In Sec.|21 we review the 
stability limits of isothermal spheres in different ensembles 
by using the turning point criterion. In Sec.|21 we consider 
specifically the grand canonical and grand microcanon- 
ical ensembles and evaluate the second order variations 
of the associated thermodynamical potential. In Sec. 0] 
we briefly discuss the connexion between thermodynam- 
ics and statistical mechanics (and field theory) for self- 
gravitating systems. In Sec. El we consider the case of gen- 
eralized thermodynamics proposed by Tsallis (1988) and 
leading to stellar polytropes (Plastino & Plastino 1997, 
Taruya & Sakagami 2002a,b, Chavanis 2002b). In Sec. 
El we discuss the stability of isothermal and polytropic 
spheres under an external pressure (Bonnor problem) . We 
provide a new and entirely analytical solution of this old 
problem. Finally, in Seed we discuss the relation between 
dynamical and thermodynamical stability of stellar sys- 
tems and gaseous spheres. 



2. Thermodynamics of self-gravitating systems 

2.1. Thermodynamical ensembles 

Consider a system of N particles, each of mass m, in- 
teracting via Newtonian gravity. Let /(r, v, i) denote the 
distribution function of the system and p{r, t) — J f(Pv 
the spatial density. The total mass and total energy can 
be expressed as 



M = / p(fr, 



(1) 



(2) 



where K is the kinetic energy and W the potential energy. 
The gravitational potential <& is related to the spatial den- 
sity by the Poisson equation 



A* = AnGp. 



(3) 



The Boltzmann entropy is given by the standard formula 
(within an additional constant term) 



S 



f 



(4) 



which can be obtained by a combinatorial analysis 
(Ogorodnikov 1965). For the moment, wc shall assume 
that the system is confined within a spherical box of ra- 
dius R so that its volume is fixed. 

In the microcanonical ensemble, the system is isolated 
and conserves energy E and particle number N (or mass 
M — Nm). Statistical equilibrium corresponds to the 
state that maximizes the entropy S at fixed E and A^^. 
Introducing Lagrange multipliers 1/T and —p./T for each 
constraint, the first order variations of entropy satisfy the 
relation 



5S 



T T 



(5) 



where T is the temperature and p the chemical potential. 
Eq. (jSJ can be regarded as the first principle of thermo- 
dynamics. 

In the canonical ensemble, the temperature and the 
particle number are fixed allowing the energy to fluctu- 
ate. In that case, the relevant thermodynamical parame- 
ter is the free energy (more precisely the Massieu func- 
tion) J = S — ^E which is related to S* by a Legendre 
transformation. According to Eq. Q, the first variations 
of J satisfy 5 J — —E5{^) — ^5N. Therefore, at statistical 
equilibrium, the system is in the state that maximizes the 
free energy J at fixed particle number N and temperature 
T. 

In the grand canonical ensemble, both temperature 
and chemical potential are fixed, allowing the energy and 
the particle number to fiuctuate. The relevant thermody- 
namical potential is now the grand potential G = S—-^E-\- 
^N. Its first variations satisfy 5G = ~E5{^) + NS{^). 
At statistical equilibrium, the system is in the state that 
maximizes G at fixed T and n/T. 
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The microcanonical, canonical and grand canonical 
ensembles are the most popular. However, we are free 
to define other thermodynamical ensembles. For exam- 
ple, Lecar & Katz (1981) have introduced a grand mi- 
crocanonical ensemble (GMCE) in which the fugacity and 
the energy are fixed. This corresponds to a thermodynam- 
ical potential JC = S + whose first variations satisfy 
6JC = ^SE+NS{^). At statistical equilibrium, the system 
is in the state that maximizes /C at fixed E and /i/T. 



2.2. The isothermal distribution 

For extensive systems, it is well-known that the statisti- 
cal ensembles described previously are all equivalent at 
the thermodynamic limit. Therefore, the choice of a par- 
ticular ensemble (in general the grand canonical one) is 
only dictated by reasons of convenience. This is not the 
case for systems with long-range interactions such as grav- 
itational systems. Indeed, the statistical ensembles are not 
interchangeable and the stability limits differ from one en- 
semble to the other. Therefore, the choice of the ensemble 
is imposed by the physical properties of the system under 
consideration. If we want to describe in terms of statis- 
tical mechanics globular clusters (Lynden-Bell & Wood 
1968) and elliptical galaxies (Lynden-Bell 1967, Hjorth 
& Madsen 1993, Chavanis & Sommeria 1998), which can 
be treated as isolated systems in a first approximation, 
the good choice is the microcanonical ensemble. If now 
the system is in contact with a radiation background 
imposing its temperature, the relevant ensemble is the 
canonical one. This ensemble may be appropriate to star 
or galaxy formation (Penston 1969, Chavanis 2002a, b), 
white dwarf and neutron stars (Hertel & Thirring 1971, 
Chavanis 2002c, d) and dark matter made of massive neu- 
trinos (Bihc & VioUier 1997, Chavanis 2002g). This canon- 
ical description is also exact in a model of self-gravitating 
Brownian particles (Chavanis et al. 2002, Sire & Chavanis 
2002). Finally, the grand canonical ensemble may be ap- 
propriate to the interstellar medium (or possibly the large- 
scale structures of the universe) , assuming that a statisti- 
cal description is relevant (de Vega et al. 1998). Possible 
applications of the grand microcanonical ensemble are dis- 
cussed by Lecar & Katz (1981). 

If we just cancel the first order variations of the ther- 
modynamical potential (under constraints appropriate to 
the ensemble considered), it is straightforward to check 
that each ensemble yields the Maxwell-Boltzmann distri- 
bution 



(6) 



Therefore, the statistical equilibrium state of a self- 
gravitating system is isothermal whatever the ensemble 
considered. The differences will come from the second 
order variations of the thermodynamical potential. The 
choice of the statistical ensemble will affect the stability of 
the system. 



The density field associated with the distribution func- 
tion JH)) is given by 



27r 



3/2 



(7) 



where we have defined the fugacity z and the inverse tem- 
perature P by the relations 



m 



(8) 



The distribution function © can then be written in terms 
of p as 



(9) 



This distribution function is a global maximum of the ther- 
modynamical potential under the usual constraints and 
for a given density field p(r) (see Padmanabhan 1989, 
Chavanis 2002a). Substituting this optimal distribution 
in Eqs. |(2J) and we get 



E = ^NkT + i / p^d^r, 



S 5N 3N ^ f2TTkT 

- = \ In 

k 2 2 



/27rfcT\ 



^In^d^r. 



(10) 



(11) 



It will be more convenient in the sequel to study the sta- 
bility problem directly from these equations, without loss 
of generality. 

2.3. The Emden equation and the Milne variables 

Inserting the relation Q in Eq. (j2Jl , we find that the gravi- 
tational potential $ is a solution of the Boltzmann-Poisson 
equation 



A$ = 47rGAe-''*. 



(12) 



If we introduce the function i}) — /?($ — $0) where <I>o is 
the gravitational potential at r = 0, the density field can 
be written 



P = Poe 



(13) 



where po is the central density. Introducing the normal- 
ized distance ^ = {AirGPpoy^^r, the Boltzmann-Poisson 
equation can be written in the standard Emden form 
(Chandrasekhar 1942) 



1 d 



(14) 



with ■i/; = '0' = Oat^ = O. This equation has to be solved 
between ^ = and S, — a corresponding to the normalized 
box radius 



a = {AnCPpof/'^R. 



(15) 
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We also recall the definition of the Milne variables that 
will be used in the sequel 



V = i^'. 



These variables satisfy the identities 

Idu I 
u at, t 



1 dv 
vd^ 



(16) 



(17) 



(18) 



which can be directly derived from Eq. (|14|l . 

The phase portrait of classical isothermal spheres is 
represented in Fig. ^ and forms a spiral. The spiral starts 
at {u, v) = (3, 0) for ^ = and tends to the limit point 
(us, Vs) = (1, 2) as ^ ^ +00. This limit point corresponds 
to the singular solution e'"^" — Ij^- For finite isother- 
mal spheres, we must consider only the part of the curve 
between ■J = and — a. 

2.4. The equilibrium phase diagram 

The thermodynamical parameters can be expressed in 
terms of the values of the Milne variables at the nor- 
malized box radius a. Writing mq = u{a) and wq = 
v{a)^ the energy and the temperature can be written as 
(Padmanabhan 1989, Chavanis 2002a): 



A 



ER 

^ GAP 

I3GM _ 
R ~ 



1 



uo 



Vq. 



(19) 



(20) 



We need also to express a in terms of the fugacity z. 
Using the expression of the gravitational potential at the 
box radius $(i?) = —GM/R and comparing Eqs. O and 
(O, we get 



3/2 



ze' 



(21) 



Expressing po in terms of a by Eq. (|15|l and introducing 
the Milne variables H16|) . we obtain 



X 



A{2TTfG^R'^z^ 



(22) 



We can also normalize the fugacity by the energy instead 
of the temperature. Combining Eqs. l|221), (EOJ and ifT^ . 

we get 



-16(27r)i°G5i?^z4£; 



-ivo 



(23) 



It is easy to show that statistical equilibrium states 
only exist for sufficiently small values of the control pa- 
rameters A, T], X and v defined previously. Above a critical 
value, there is no possible equilibrium state and the system 
is expected to collapse. For example, for A > Ac = 0.335 




Fig. 1. This graphical construction, first introduced by 
Padmanabhan (1989) in MCE, shows very simply the ex- 
istence of an upper bound for A, r\ and x above which there 
is no hydrostatic equilibrium for an isothermal gas (dotted 
line: A > Ac; solid line: A = Ac; dashed line: A < Ac with 
similar convention for 77 and x). 



(which corresponds to sufficiently negative energies), we 
have the well-know gravothermal catastrophe (Antonov 
1962, Lynden-Bell & Wood 1968). Similarly, for tj > r]c = 
2.52 (which corresponds to small temperatures or large 
masses), we have an isothermal collapse (Chavanis 2002a, 
Chavanis et al. 2002). We can show the existence of such 
bounds by a simple graphical construction, without nu- 
merical work. Indeed, the equilibrium structure of the sys- 
tem (characterized by its degree of central concentration 
a) is determined by the intersection between the spiral in 
the {u, v) plane and the curves defined by Eqs. ifT^ . 
and (|22|l . If there is no intersection, the system cannot be 
in statistical equilibrium. This graphical construction is 
made exphcitly in Fig. [H in MCE, CE and GCE. 

In the microcanonical ensemble, the control parame- 
ters are the energy E and the particle number N. The 
parameter conjugate to the energy (with respect to the 
entropy S) is the inverse temperature /3. The curve giving 
the inverse temperature as a function of energy for a fixed 
particle number N is drawn in Fig. [3 The parameter con- 
jugate to the particle number N is the ratio —fi/T, related 
to the fugacity z. In order to represent the fugacity as a 
function of A^, for a fixed energy, we use Eqs. H19() and (|23() . 
The —p./T vs N curve can then be deduced from Figs.|21 
and 01 According to standard turning point arguments 
(Katz 1978), the series of equilibria becomes unstable at 
the point of minimum energy (for a given mass) or at the 
point of minimum mass (for a given negative energy). The 
condition dA/da = implies (Padmanabhan 1989) 



4mq -I- 2uoVo — llup -f 3 = 0. 



(24) 
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-0.4 




A=-ER/GM 



Fig. 2. Temperature vs energy plot at fixed particle num- 
ber. The series of equilibria becomes unstable at MCE in 
the microcanonical ensemble and at CE in the canonical 
ensemble. 
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Fig. 3. Fugacity vs particle number at fixed energy. The 
series of equilibria becomes unstable at GMCEl in the 
grand microcanonical ensemble. 



The values of a at which the parameter A is extremum are 
determined by the intersections between the spiral in the 
{u, v) plane and the parabole defined by Eq. (|24|l . see Fig. 
O The change of stability occurs for the smallest value 
of a. Its numerical value is a = 34.4 corresponding to a 
density contrast of 709. Gravitational instability occurs in 
MCE when the specific heats Cy = {dE / dT)M.v becomes 
zero, passing from negative to positive values. 

In the canonical ensemble, the control parameters are 
the inverse temperature (3 and the particle number N . 
The parameter conjugate to the inverse temperature (with 
respect to the free energy J) is —E. The curve giving 
—E as a, function of [3 for a fixed particle number can 
be deduced from Fig. |21 The parameter conjugate to the 
particle number is — ^/T, related to the fugacity z. We 
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Fig. 4. Enlargement of Fig.Onear the spiral. In the grand 
microcanonical ensemble, a change of stability occurs at 
GMCE2 and GMCE3. Between these points the series of 
equilibria is stable again. After GMCE3, the series be- 
comes and remains unstable. 
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Fig. 5. Graphical construction determining the turning 
points of the control parameter at which a change of sta- 
bihty occurs in MCE, CE and GCE. 

can represent the fugacity as a function of mass, for a 
given temperature, by using the relations (|22() and (|20(l . 
The corresponding curve can be deduced from the ry vs 
X plot represented in Fig. El In the canonical ensemble, 
the series of equilibria becomes unstable at the point of 
minimum temperature (for a given mass) or at the point of 
maximum mass (for a given temperature). The condition 
dr] / da = is equivalent to (Chavanis 2002a) 



uo 



1. 



(25) 



The values of a at which the parameter 77 is extremum are 
determined by the intersections between the spiral in the 
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Fig. 6. Fugacity vs mass at fixed temperature. The se- 
ries of equilibria becomes unstable at GCE in the grand 
canonical ensemble. 
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Fig. 7. Energy vs temperature at fixed fugacity. 

{u,v) plane and the straight line defined by Eq. (|25|l . see 
Fig-El The change of stability occurs for the smallest value 
of a. Its numerical value is a = 8.99 corresponding to a 
density contrast of 32.1. Gravitational instability in CE 
occurs when the specific heats becomes negative passing 
by Cy = oo. 

In the grand canonical ensemble, the control param- 
eters are the inverse temperature (3 and the ratio fi/T 
related to the fugacity z. The parameter conjugate to the 
inverse temperature (with respect to the grand potential 
G) is —E. The curve giving — iJ as a function of f3 for a 
fixed fugacity z is determined by Eqs. and and 
can be deduced from Figs. [7| and |S1 The parameter conju- 
gate to the ratio /i/T is the particle number N. The curve 
giving TV as a function of z for a given temperature can 
be deduced from Fig. In the grand canonical ensem- 
ble, the series of equilibria becomes unstable at the point 
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Fig. 8. Enlargement of Fig.[7|near the spiral. 
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Fig. 9. Evolution of x as a function of the central concen- 
tration a. In dashed line we have represented the asymp- 
totic value Xs — 4e~'* corresponding to the singular solu- 
tion. The A{a) and ri{a) curves have a similar behaviour. 

of maximum temperature (for a given fugacity) or at the 
point of maximum fugacity (for a given temperature) . The 
condition dx/da — implies 



uqVo = 2. 



(26) 



The values of a at which the parameter x is extrcmum 
(see Fig. ^ are determined by the intersections between 
the spiral in the (u, v) plane and the hyperbole defined 
by Eq. (|25|l . see Fig.|Sl The change of stability occurs for 
the smallest value of a. Its numerical value is a — 1.78 
corresponding to a density contrast of 1.58. 

In the grand microcanonical ensemble, the control pa- 
rameters are the energy E and the ratio /i/T related to the 
fugacity z. The parameter conjugate to the energy (with 
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Fig. 10. Evolution of as a function of the central concen- 
tration a. In dashed line we have represented the asymp- 
totic value Vg = 2^e^^ corresponding to the singular solu- 
tion. Contrary to Fig.|51 the second extremum a2 is above 
the asymptote v — Vs- This appears to be the reason for 
the regain of stability between a2 and a^. 



respect to the potential /C) is /3. The curve giving /3 as a 
function of E for a fixed fugacity z can be deduced from 
Figs. Hand |S| The parameter conjugate to the ratio fi/T 
is the particle number A^. The N vs fi/T plot for a given 
energy can be deduced from Figs. 01 and 0] In the grand 
microcanonical ensemble, the series of equilibria becomes 
unstable at the point of maximum energy (for a given fu- 
gacity) or at the point of maximum fugacity (for a given 
positive energy). The condition dv/da = implies 



BuqUo — IOwqWo — 17uo + 21 = 0. 



(27) 



The values of a at which the parameter v is extremum 
(see Fig. I1U|I are determined by the intersections between 
the spiral in the (u, v) plane and the curve defined by 
Eq. H27|) . see Fig.^] The series of equilibria becomes un- 
stable for ai — 1.89 (density contrast 1.66) which corre- 
sponds to the first intersection. At that point, the v vs 
X curve rotates clockwise (see Fig. Ejl. However, at the 
second turning point 02 = 5.95 (density contrast 11.6) 
the v vs X curve rotates anticlockwise (see Fig. |S} so 
that stability is regained (this is proved in Katz 1978). 
Stability is lost again at — 13.8 (density contrast 92.7) 
and never recovered afterwards. This curious behaviour in 
GMCE was first noted by Lecar & Katz (1981). A sim- 
ilar lost-gain-lost stability behaviour also occurs for self- 
gravitating fermions and hard sphere models in MCE and 
CE (Chavanis & Sommeria 1998, Chavanis 2002c). In that 
case, it is related to the occurence of phase transitions be- 
tween "gaseous" and "condensed" states. 



a,=13.8 




Fig. 11. Graphical construction determining the turning 
points of the control parameter v in GMCE. The series of 
equilibria is stable for a < ai, unstable for ai < a < 012, 
stable again for a2 < a < and unstable for a > a^. 

3. Thermodynamical stability of self-gravitating 
systems 

We now determine the condition of thermodynamical sta- 
bility by computing the second order variations of the 
thermodynamical potential explicitly. This study has al- 
ready been done in MCE (Padmanabhan 1989) and CE 
(Chavanis 2002a). In the following, we perform the anal- 
ysis in GCE and GMCE. 

3.1. The grand canonical ensemble 

The second order variations of the grand potential G are 



m 



2p 2T 



SpS^d^r. 



(28) 



An isothermal sphere is stable in the grand canonical en- 
semble if ^^G < for all variations 6p of the density pro- 
file. We shall restrict ourselves to spherically symmetric 
perturbations since only such perturbations can trigger 
gravitational instability for non rotating systems (Horwitz 
& Katz 1978). In the grand canonical ensemble, the mass 
of the system is not conserved. Therefore, the potential at 
the wall changes according to 



(5$(i?) = - 



GSM 
R 



(29) 



On the other hand, according to the Gauss theorem, we 
have 



dS^ . ^, G5M 



(30) 



These boundary conditions suggest the introduction of the 
new function 



X(r) = rS^{r), 



(31) 
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which satisfies 
X{0) = 0, X"iO) 



0, 



X'{R) = 0, 



(32) 



where we have used Eqs. (|29|) and H30|l to obtain the third 
equality (the second equality is obvious for a spherically 
symmetric system). 

The variations (5$ and 6p are related to each other by 
the Poisson equation 

r-^ dr \ ar J 

In terms of the new variable (|31|l . Eq. H33I) becomes 

— = AirGSp. (34) 
r 

Expressing the second order variations of G in terms of 
X, we get 



S'G = 



k 



ll\2 



-dr 



1 



X"Xdr. (35) 



m87rG2 7o p ' 2TG 
Integrating by parts, we can put Eq. H35() in the form 



5^G 



1 

2G2 





G 


k d 1 


^ 1 d\ 


/ drX' 








'o 




m dr ' 


\^ Airp dr ) 



X'. 



(36) 



The boundary terms arising from the integration by parts 
are seen to vanish by virtue of Eq. H32|l . We can therefore 
reduce the stability analysis of isothermal spheres in the 
grand canonical ensemble to the eigenvalue problem 



k d 1 


^ 1 d\ 


G" 


m dr ' 


yiirp dr J 


+ r 



X',{r) = XX',{r). 



(37) 



If all the eigenvalues A are negative, then the isothermal 
sphere is a maximum of the thermodynamical potential G 
(and is therefore stable in GCE). If at least one eigenvalue 
is positive, it is an unstable saddle point. The condition 
of marginal stability A = determines the transition be- 
tween stability and instability. We thus have to solve the 
equation 



k d / 1 dX' 



GX' 
T 



0. 



(38) 



m dr \ 47rp dr 

Integrating once and using the boundary conditions (|32|l 
at r = 0, the constant of integration is seen to vanish and 
we get 

GX 



m 47rp T 



= 0. 



(39) 



Introducing the dimensionless variables defined in Sec. 12. 31 
we arc led to solve the problem 



e^X" + X = 0, 
X{0) = X'{a) = 0. 
Let us denote by 



de ' 



(40) 
(41) 

(42) 



the differential operator which occurs in Eq. (|40|l . Using 
the Emden equation 114() . we readily check that 

m^t c{e^') = 2^, /:(e^) = e+^^. (43) 

Therefore, the general solution of Eq. ()40|l which satisfies 
the condition X = a.t the origin is 



XiO^c.a^yj' -2), 



(44) 



where ci is an arbitrary constant. The boundary condition 
X'{a) = implies 

Q,2g-v(") _ 2 = 0. (45) 

This relation determines the value of a at which the se- 
ries of equilibria becomes unstable in the grand canoni- 
cal ensemble. In terms of the Milne variables, Eq. H45|l is 
equivalent to 



uoVq = 2, 



(46) 



which is precisely the condition (|26() obtained with the 
turning point criterion. We are now in a position to deter- 
mine the structure of the perturbation profile that triggers 
the instability in GCE. Using Eq. H34|l . we have 



(47) 



(48) 



(49) 



5p _ X" 
po 47r^ ' 

with Eq. (gH) for X. This yields 

Sp^ 
Po 47r' 

or, introducing the Milne variables, 

^ = f(2-.). 
p An 

We note that Eq. H49|l coincides with the expression found 
in the canonical ensemble (Chavanis 2002a). The number 
of nodes in the profile 6p is determined by the number 
of intersections between the spiral in the (u, v) plane and 
the line v — 2. We see in Fig. |S1 that there is no intersec- 
tion for ^ < a = 1.78, where 1.78 is the stability limit in 
GCE. Therefore, collecting the results obtained by differ- 
ent authors, the perturbation that triggers gravitational 
instability has two nodes in MCE ( "core- halo" structure) , 
one node in CE and no node in GCE (see Fig. I12|l. 

3.2. The grand microcanonical ensemble 

The grand microcanonical ensemble is interesting in its 
own right but also because it exhibits a lost-gain-lost sta- 
bility behaviour similar to that found for self-gravitating 
fermions and hard sphere models (Chavanis 2002c). Now, 
the stability analysis can be conducted analytically in 
GMCE for classical isothermal spheres while only numer- 
ical results are available for fermions and hard sphere sys- 
tems. On a technical point of view, GMCE is the most 
complicated and richest of all thermodynamical ensem- 
bles. We can argue, however, that GMCE may not be re- 
alized in practice. It is not clear, indeed, what mechanism 
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0.5 



-0.5 







{GCE)^\ 
(MCE) 


(CE) 



The stability analysis reduces therefore to the study of the 
eigenvalue equation 

dr'K{r, r')X';^{r') = XX';^{r). (55) 

The condition of marginal stability (A = 0) reads 



GX' 
T 



k d / 1 dX' 
m dr \47r/9 dr 

X i \{<^ 



3kT 
'2m 







3kT 
2m 



0.2 



0.4 0.6 
r/R 



0.8 



Integrating once, we obtain 

GX k X" _ 2 
~T~ ^mA^p^ 3NkT^ 



X'dr. 



3kT\ 
2m J 



(56) 



Fig. 12. First mode of instability in MCE, CE and GCE. 



3kT 
2m 



X'dr. 



(57) 



The gravitational potential <i> is related to the Emden 
function by the relation ■(/; = /?($ — <I>o). Using 

can be devised so that two systems can share particles <i>(i?) = ~GM/R and Eq. j^OJl, we get 



without also sharing energy. One example might be the 
case of relativistic particles created and destroyed at equi- 
librium in a cosmological settling. 

In GMCE, we have to maximize the thermodynamical mensionless variables defined in Sec.E21 we find that 
potential 



/3$ = -0 - V'(a) - vq. (58) 
Inserting this relation in Eq. 157|l and introducing the di- 



5 3, /2TTkT 

- + - In 

2 2 V ™ 

^In^dV 



/C 
k 



m m 
at fixed energy 
3 



kT 



X + e'I'X" = 



E 



-kT 



P-d'v- 



p^d^r. 



(50) 



(51) 



1 
2 



— ^^:{a) — vo + 



X^e-'^d^ + Xia) }, 







(59) 



where we have integrated by parts and used the identity 
(^■0)" = ^e~'^ equivalent to the Emden equation (|14|) . 
Writing 



V = 



pa 

X{a) - J X^e-'^d^ 



We vary the density profile and use Eq. (|51|l to determine 
the corresponding variation of temperature. After some 
algebra, we find that the second order variations of the 
thermodynamical potential can be written 



W 



V, 



S^JC = - — 
2T 

1 



SpS^d'^r 

m 



3NkT^ 



3kT\ 
2m 



M 

2p 



jSp d^ 



■d^r 



Introducing the variable X defined previously and inte- 
grating by parts, we can put the problem in the form 



with 



1 

2G2 



K{r, r') 



dr / dr'X'{r)K{r,r')X'{r'), 



we are led to solve the problem 
e'^X" + X = Vipi + T^f , 
X{Q) = X'{a) = 0. 
(52) Seeking a solution of Eq. H()2|) in the form 
X ^ai£_+ a2CV + 03^-0, 
and using the identities (|43|l . we find that 
X = {W-V)S,+ bif^' - 20 + V^i;. 



(53) 



(60) 
(61) 

(62) 
(63) 

(64) 

(65) 



The constant b is determined by the boundary condition 
X'{a) — 0. We thus find that the general solution of Eqs. 
(|?^ and jSnj is 



3NkT^ 



G k d f 1 d 
T m dr \4ttp dr 
3kT 



2m 



(r) 



d{r — r') 
3kT\ 



2m 



r 



(54) 



■0(a) 
V 



Vq - 
1 



2 uoVq — 2 



(66) 
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where we have used Eq. (|61l) . The values of a at which a 
change of stabiUty occurs (A = 0) are obtained by inserting 
the solution Ht)t)|) for X in Eq. (|60|l . We need the identities 



^^^e-'^d^ = avoWa) + vq - 6 + 2uo). 



(67) 
(68) 
(69) 



which are derived in Appendix 1X1 After simplification, we 
obtain the condition 



IOwqWo - 17mo + 21 = 0, 



(70) 



which coincides with the condition (|27|l deduced from the 
turning point criterion. The perturbation profile at the 
critical points at which A = is given by Eq. H47(l with 
Eq. ^ for X. Using Eqs. (O, (E3) and JMI), we get 



— = -7!-(2uoWo -6 + v), 
p 47r 



(71) 



where ci is a constant. The nodes of the perturbation pro- 
file are determined by the condition 



(72) 



with < a. The number of nodes can be found by a 
graphical construction. In Fig. E| we plot the lines de- 
fined by Eq. 172|) in the (w, v) plane for the four first crit- 
ical points in the series of equilibria. The corresponding 
perturbation density profiles are plotted in Fig. ^1 The 
first mode of instability ai has no node. The profile of the 
second mode a2 increases (in absolute value) with dis- 
tance. The absence of node at a2 is a signature of the fact 
that stability is regained at that point. The mode as at 
which stability is lost again has one node, the mode 
two nodes etc... 




Fig. 13. Graphical construction determining the number 
of nodes for the different modes of instability in GMCE. 
The horizontal lines correspond to the condition H72(l for 
the four first critical points. 




4. Connexion with statistical mechanics 

In this section, we briefly discuss the relation between the 
thermodynamics and the statistical mechanics (and field 
theory) of self-gravitating systems. We shall just discuss 
the physical ideas without entering into technical details. 

In the microcanonical ensemble, the object of funda- 
mental interest is the density of states 

N 

g{E) = 5{E-H)Y[ d'rd^v,, (73) 
where 

is the Hamiltonian of the self-gravitating gas. By defini- 
tion, g{E)dE is proportional to the number of microstates 



Fig. 14. Density perturbation profiles for the four first 
critical points in the series of equilibria in GMCE. We 
have fixed the constant ci such that Sp/ p = 1 at r = 0. 

(specified by the position and the velocity of all N parti- 
cles) with energy between E and E + dE. The entropy is 
defined by 

S{E)^\ng{E). (75) 

As is customary in statistical mechanics, we decompose 
phase space into macrocells and divide these macrocells 
into a large number of microcells. A macrostate is char- 
acterized by the number of particles in each macrocell 
(irrespective of their precise position and velocity in the 
cell), or equivalently by the smooth distribution function 
/(r,v). If we call Wi the number of microstates corre- 
sponding to the macrostate i, we can rewrite Eq. I|73(l 
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formally in the form 

giE)= J2W.^Y.^''-I: ^"''^ (76) 

Ei=E Ei = E Ei=E 

where the sum runs over all macrostates with energy E 
(and mass M). In addition, Si — hiWi is the Boltzmann 
entropy I0J of the macrostate i (i.e., 5*^ = S[f]) and 
Si = Si/N the entropy per particle. The thermodynamic 
limit for self-gravitating systems corresponds to — > +oo 
in such a way that all the dimensionless control parame- 
ters defined in Sec. 12.41 remain finite. This can be written 
for example N,R^ +00 with N/R finite (and E/N r^l, 
/3 ~ 1). This is a very unusual thermodynamic limit due 
to the long-range nature of the interactions. Now, at the 
thermodynamic limit, the sum in Eq. H76(l is dominated 
by the macrostate /* which maximizes the Boltzmann en- 
tropy 101) at fixed mass and energy: 

g{E)^e^^f'\ S{E)^S[U]. (77) 

It has to be noted that the sum in Eq. (|76(l is sharply 
peaked around which corresponds to what is called a 
large deviation property in mathematics. Therefore, the 
mean-field approximation is exact for self-gravitating sys- 
tems at the thermodynamic limit. 

In fact, if we look at the problem more carefully, the 
argument leading to Eq. H77|) is not so obvious because 
the entropy S[f] has no global maximum, except if we in- 
troduce small-scale and large-scale cut-offs. The introduc- 
tion of these cut-offs is necessary to define correctly the 
density of state (|73|l . otherwise the integral would diverge 
(Padmanabhan 1990). The divergence at large scales is as- 
sociated to the fact that a stellar system has the tendency 
to evaporate under the effect of encounters. In reality, stel- 
lar systems are not allowed to extend to infinity because 
they interact with the surrounding. For example, globular 
clusters are subject to the tides of a nearby galaxy so that 
the Boltzmann distribution has to be truncated at high 
energies. The Michie-King model, 

r ^(e-/5^_e-/3^") 


which is a truncated isothermal, can take into account 
the evaporation of high energy stars and provides a good 
description of about 80% of globular clusters (Binney 
& Tremaine 1987). It can be derived from the Fokker- 
Planck equation (taking into account the encounters be- 
tween stars) by imposing that the distribution function 
vanishes at the escape energy e = e^. In that case, the 
system is not truly static since it gradually loses stars 
but we can consider that a globular cluster passes by a 
succession of quasi-equilibrium configurations. Instead of 
working with truncated models, it is more convenient to 
confine the system within a box of radius R, the box radius 
playing the role of a tidal radius. The divergence at small 
scales of the density of states reflects the natural tendency 
of a stellar system to form binaries. In this extreme case, 
the size a of the stars provides a small-scale cut-off. The 



influence of a small-scale cut-off on the form of the equi- 
librium phase diagram is discussed by Chavanis (2002c) 
and Chavanis & Ispolatov (2002). When the small-scale 
cut-off a — !■ 0, which is the limit relevant for globular clus- 
ters since the size of the stars is in general much smaller 
than the interstellar distance, the entropy S[f] possesses 
two maxima (for sufficiently high energies). The local en- 
tropy maximum describes a "gaseous" phase insensitive 
to the small scale cut-off. It corresponds to the isother- 
mal spheres studied previously. The global entropy max- 
imum is dominated by the small scale cut-off and corre- 
sponds to a "condensed" phase made of a few stars close 
together (or simply a single binary) surrounded by a dif- 
fuse halo. This "core-halo" structure is the most probable 
configuration of a stellar system so it is expected to be 
reached at equilibrium. Since the observed globular clus- 
ters, described by the Michie-King model, do not possess 
this core-halo structure, we conclude that they are not in 
true equilibrium but that they are slowly evolving in that 
direction. However, the formation of binaries can take ex- 
tremely long times, much larger than the age of the uni- 
verse. Indeed, equilibrium statistical mechanics tells noth- 
ing about the relaxation time. For the timescales contem- 
plated in astrophysics, globular clusters can be considered 
as metastable equilibrium states corresponding to local en- 
tropy maxima. Indeed, the probability of transition from 
a gaseous state to a condensed state (corresponding to 
the formation of a binary) is extremely small (Katz & 
Okamoto 2000, Chavanis & Ispolatov 2002). Therefore, in 
the calculation of the density of states H73|l . we must dis- 
card by hands the states which correspond to collapsed 
configurations (binaries) since they cannot be reached in 
the timescales of interest. Accordingly, Eq. I|77|) is correct 
with being the local maximum of S[f]. This is precisely 
the situation studied in Secs.ElandlHl In fact, this picture 
is correct only for globular clusters with high energy E, 
well above the Antonov limit Ec- Because of evaporation, 
the energy of a globular cluster slowly decreases until it 
reaches the point at which the system cannot be in equilib- 
rium anymore (e.g., Katz 1980). At that point, an isother- 
mal sphere becomes unstable and the system undergoes 
a gravitational collapse (gravothermal catastrophe). The 
evolution proceeds self-similarly with the formation of a 
high density core with a shrinking radius and a small mass 
(Henon 1961, Cohn 1980, Lynden-BeU & Eggleton 1980). 
This core collapse concerns typically 20% of globular clus- 
ters. In theory, the central density becomes infinite in a 
finite time. In practice, the formation of binaries can re- 
lease sufficient energy to stop the collapse (Henon 1961) 
and even drive a reexpansion of the system (Inagaki & 
Lynden-BeU 1983). Then, a series of gravothermal oscilla- 
tions should follow (Bettwieser & Sugimoto 1984). 

In the canonical ensemble, the object of interest is the 
partition fmiction 

N 

Z{P)= / e-'^^n^'^'^'^- (79) 
■' i=i 
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The free energy is defined by 



1 



F{[3)=--\nZ. 



(80) 



The partition function is related to the density of states 
by the Laplace transform 



/i-oo 
dE g{E)e-^^. 
-OC 

Using Eq. (jTGfl . it can be rewritten 



.Si-fJEi 



.Ji 



(81) 



(82) 



where Ji = Si — fiEi is the Massieu function and the 
sum runs over all macrostates with mass M. In the ther- 
modynamic limit, the partition function is dominated by 
the contribution of the macrostate which maximizes the 
Massieu function Ji = J[f] at fixed mass. In this limit, 
the mean field approximation is exact and we have 



Z(/3) 



F = E -TS. 



(83) 



Again, the Massieu function J[f \ has no global maximum 
in the absence of cut-off. If we introduce a small scale cut- 
off a and a large scale cut-off R and consider the limit 
a/R <C 1, it is found that the Massieu function J[/] has 
two maxima (for sufficiently high temperatures). The local 
maximum corresponds to isothermal "gaseous" configura- 
tions and the global maximum to a "condensed" structure 
in which all the particles are packed together (Aronson 
& Hansen 1972, Chavanis 2002c, Chavanis & Ispolatov 
2002). As a ^ 0, the density profile of the condensed 
state tends to a Dirac peak (Kiessling 1989). As in the 
microcanonical ensemble, the "gaseous" configurations are 
metastable but they can be long lived and relevant for as- 
trophysical purposes. This justifies why we only select the 
contribution of the local maximum of J[f] in Eq. H83(l . 
This treatment is, however, valid only for sufficiently high 
temperatures. Below T^, a phase transition similar to the 
gravothermal catastrophe occurs in the canonical ensem- 
ble and leads to a condensed structure containing almost 
all the mass. The inequivalence of statistical ensembles re- 
garding the formation of binaries (in MCE) or Dirac peak 
(in CE) is further discussed in Appendices A and B of Sire 
& Chavanis (2002). 

Finally, the grand canonical partition function is de- 
fined by 



Zgc — 



Af=0 



N r ^ 

i—1 



(84) 



where z is the fugacity. Applying the standard Hubbard- 
Stratanovich transformation, it can be rewritten in 
the form of a path integral (Horwitz & Katz 1978, 
Padmanabhan 1991, de Vega & Sanchez 2002) 



Zgc 



V(t> e 



(85) 



47r- 



Grr? 



T 



'.GttJ/'^Vt. 



for the Liouville action 



(86) 



(87) 



It should be emphasized that (/) is a formal field which is 
not the gravitational field. In the mean-field approxima- 
tion, the path integral is dominated by the contribution 
of the field 4>o which maximizes the Liouville action (|87|l . 
The cancellation of the first order variations gives 

A(/.o = -//V°, (88) 

which is similar to the Boltzmann-Poisson equation (|12|) . 
This allows one to identify 0o with the equilibrium gravi- 
tational potential $ (up to a negative proportionality fac- 
tor). The second order variations of A at the critical point 
must be negative for A[0o] to be a maximmn. This yields 
the condition 



J d^r S(j){A + n^e'l"'}S(f><0, 



(89) 



which is equivalent to the maximization of the grand po- 
tential G (see Sec. I3.1|l . Therefore, the thermodynamical 
approach is equivalent to the statistical mechanics or field 
theory approach. Wc note, however, that the Liouville ac- 
tion A[(j)\ is not the same functional as the grand potential 
G[p] although they give the same critical points and the 
same conditions of stability. 

5. Generalized thermodynamics and Tsallis 
entropy 

5.1. General considerations 

It has been recently argued that the classical Boltzmann 
entropy may not be relevant for non-extensive systems 
and that Tsallis entropies, also called g-entropies, should 
be used instead. These entropies can be written 



So 



(90) 



where g is a real number. For q —^ 1, this expression re- 
turns the classical Boltzmann entropy as a particular 
case. 

The first application of Tsallis generalized thermody- 
namics in stellar systems was made by Plastino & Plastino 
(1997) who showed that the extremization of Tsallis en- 
tropy leads to stellar polytropes with an index 



3 

n = - + - 
2 q 



1 



1 



(91) 



For q — > 1, or n — > -|-oo, we recover isothermal spheres. 
For q > 9/7 (i.e. n < 5), the density drops to zero at a 
finite distance so that the mass is finite. Then, Taruya & 
Sakagami (2002a) considered the stability problem in the 
microcanonical ensemble by extending Padmanabhan's 
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classical analysis of the Antonov instability to the case of 
polytropic distribution functions. A configuration is sta- 
ble in the sense of generalized thermodynamics if it cor- 
responds to a maximum of Tsallis entropy at fixed mass 
and energy. They considered polytropic spheres confined 
within a box and showed that polytropes with index n > 5 
become unstable above a certain density contrast while 
polytropes with n < 5 are always stable. On the other 
hand, if we consider the canonical ensemble, we have to 
select maxima of Tsallis free energy Jq = Sq — (3E at fixed 
mass and /? (the generalized inverse temperature). It is 
found (Taruya & Sakagami 2002b; see also Sec. 15. 5|) that 
polytropes with index rt > 3 become unstable above a 
certain density contrast while polytropes with n < 3 are 
always stable. In our sense, the main interest of the q- 
entropies (and their attractive nature) is to offer a simple 
generalization of classical thermodynamics which leads to 
models that are still analytically tractable since exponen- 
tials are replaced by power laws. This approach is interest- 
ing to develop, at least formally, as it offers a nice connex- 
ion between polytropic and isothermal spheres, which are 
the most popular mathematical models of self-gravitating 
systems. We shall therefore devote a section to this gen- 
eralization. However, in Seed we shall argue that Tsallis 
entropies are just a particular case of i/-functions S[f] 
whose maximization at fixed mass and energy determines 
nonlinearly stable stationary solutions of the Vlasov equa- 
tion (Tremaine et al 1986). We shall also argue that the 
maximization of the J- function J[f] — S — (3E at fixed 
mass determines nonlinearly stable stationary solutions 
of the Euler- Jeans equations. This discussion should de- 
mystify the concept of generalized thermodynamics intro- 
duced by Tsallis (1988). 

5.2. Stellar polytropes 

The extremization of Tsallis entropy at fixed mass and en- 
ergy yields (Plastino & Plastino 1997, Taruya & Sakagami 
2002a) 



$ 

2 



A = 



(<?-!)/? 



1 



l + {q-\){^i/T) 



(92) 



(93) 



where f3 — 1/T and fi/T (related to a) are Lagrange 
multipliers which can be called inverse temperature and 
chemical potential in the generalized sense. The distribu- 
tion function (|92|l corresponds to stellar polytropes that 
were first introduced by Plummer (1911). Their derivation 
from a variational principle was discussed by Ipser (1974) 
among others. Note that Eq. is vahd for < 2(a— $). 
For > 2(q; — <&), we set / = 0. The spatial density 
p = J f(P'v and the pressure p — ^ J fv^(Pv can be ex- 
pressed as 



p = 4V27rA(a - $)"B(3/2, n - 1/2), 



(94) 



AV27TA{a - $; 



ri+l 



n + 1 



5(3/2, n- 1/2), 



(95) 



with B{a,b) being the (3 function and use has been made 
of Eq. (|91() . Furthermore, in obtaining Eq. (|95|) . we have 
used the identity B{m + l,n) — mB{m, n) / [ni + n). The 
integrability condition J /d'^v < oo requires that n > 1/2. 
Eliminating the gravitational potential between the rela- 
tions (|94|) and H95|l . we recover the well-known fact that 
stellar polytropes satisfy the equation of state 



p — Kp'' 



7=1 



1 

n 



(96) 



like gaseous polytropes (note, however, that these systems 
do not have the same phase space distribution function, 
unlike isothermal spheres). In the present context, the 
polytropic constant is given by 



K 



1 



(n + 1) 



4%/27rAB(3/2,n- 1/2) 



-l/n 



(97) 



Using Eqs. H92() and H94() . the distribution function can 
be written as a function of the density as 



with 



l/n 



i;V2 



{n + l)K 



n-3/2 



Z = 4V27rB(3/2,n- l/2)[ii'(n+ 1)] 



3/2 



(98) 



(99) 



Using the foregoing relations, it is possible to express 
the total energy and the entropy in terms of p and p ac- 
cording to (Taruya & Sakagami 2002b; see also Sec. I7.8|l 



S 



pd^r - M 



(100) 



(101) 



5.3. The Lane-Emden equation 

The equilibrium structure of a polytrope is determined by 
substituting the relation (|94|l between p and the gravita- 
tional potential $ inside the Poisson equation 0. This 
is equivalent to using the condition of hydrostatic equilib- 
rium Wp ~ — /3V$ with the equation of state 196|) . Letting 



P = Pot 



K{n+l) 



1/2 



(102) 



where po is the central density, we can reduce the condition 
of hydrostatic equilibrium to the Lane-Emden equation 
(Chandrasekhar 1942) 



,d9_ 
di 



(103) 



with boundary conditions d{0) = 1, ^'(0) = 0. 

For n < 5, the density vanishes at some radius iden- 
tified with the radius of the polytrope (in terms of the 
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normalized distance, this corresponds to = at ^ = ^i). 
Such polytropic configurations (with infinite density con- 
trasts) will be refered to as complete polytropes. We shall 
also consider the case of incomplete polytropes, with arbi- 
trary index, enclosed within a box of radius R. For these 
configurations, the Lane-Emden equation must be solved 
between ^ = and ^ = a with 



-l/n 



K{n + 1) 



1/2 



R. 



(104) 



We also recall the expression of the Milne variables 
that will be used in the sequel. 



^0' 



These variables satisfy the identities 

Idu \ 
u at, 5 



Idv 1 , 



1), 



(105) 



(106) 



(107) 



which can be derived from Eq. (|103|1 . The phase portrait 
of polytropes in the (u, v) plane and the properties that 
we shall need in the following are summarized in Chavanis 
(2002b). 



^ - 




Fig. 15. Fugacity vs mass (for a given temperature) in 
the framework of generalized thermodynamics for different 
polytropic index. For 1 < n. < 3, x has one turning point 
but not rj] for 3 < n < 5, both % and rj have one turning 
point; for n > 5, % and rj have an infinity of turning points. 
For n < 5, complete polytropes correspond to the terminal 
point of the rj — x curve. For n > 5, the curve has a spiral 
behaviour towards the limit point ('7s, Xs) corresponding 
to the singular sphere [us^Vs] = (^, ;^). 



5.4. Turning point criterion 

Eliminating A between Eqs. (193(1 and 1)9 7() . we find that 



where we have defined (Taruya & Sakagami 2002b) 
(n-i)"- 



Q = 



16\/27r2B(3/2,n- l/2)(7- 



1) 



(108) 



(109) 



In a preceding paper (Chavanis 2002b), we have estab- 
lished that the mass of the configuration is related to the 
central density po (through the parameter a) by the rela- 
tion 



M 



AttG 



K{1 



1 



R- 



(110) 



We had conjectured that -q plays the role of a tempera- 
ture in the context of Tsallis generalized thermodynam- 
ics. This is indeed the case since, using Eq. ()108|l . rj can 
be expressed in terms of [3 as 



1 f G"M"-i 



{n + 1)" 



(111) 



For n — > +O0, this parameter becomes equivalent to the 
parameter r] = (3GM/R introduced in the study of isother- 
mal spheres (Chavanis 2002a). 

In the framework of generalized thermodynamics, the 
parameter conjugate to the energy with respect to Tsallis 



entropy Sq is the inverse temperature /3 (microcanoni- 
cal description) and the parameter conjugate to the in- 
verse temperature with respect to Tsallis free energy 
Jq = Sq — j3E is —E (canonical description). The E ~ T 
curve has been plotted and discussed in Chavanis (2002b) 
for different values of the polytropic index. According to 
the turning point criterion, the series of equilibria becomes 
unstable in the canonical ensemble for drj/da = 0. We have 
shown that this condition is equivalent to 



ri — 3 
n — 1 



(112) 



For n < 3, the polytropes are always stable since the 77(a) 
curve is monotonic. The polytropes start to be unstable 
(for sufficiently high density contrasts) when rt > 3. This 
generalized thermodynamical stability criterion coincides 
with the Jeans dynamical stability criterion based on the 
Navier-Stokes equations for a polytropic gas (Chavanis 
2002b). Since K is considered as a constant when we an- 
alyze the stability of polytropic spheres with respect to 
the Navier-Stokes equations, and since K is related to the 
temperature via Eq. H108|) . this corresponds to a sort of 
"canonical situation" . This explains why the condition of 
dynamical stability for gaseous polytropes is equivalent 
to the condition of thermodynamical stability for stellar 
polytropes in the canonical ensemble. This equivalence is 
further discussed in Sec. H 
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We now turn to the grand canonical situation. 
Combining Eqs. ((M)), and H109|) . we can write the 
density field as 



Q 



l-n 



47r(n+ l)"/33/2 



T 



/3$ 



(113) 



This expression suggests to define a generalized fugacity 
by the relation n z^l"^ = n — 3/2 + /i/T so that z is equiv- 
alent to the usual fugacity e~ when n +00. Taking 
r = i? in Eq. (|113|l and comparing with Eq. I|1Q2|I . we get 



Att{ti + l)"/33/2 



R ) ' 



(114) 



where we have used ^{R) — —GM / R. Expressing the cen- 
tral density in terms of a, using Eqs. H104|l and (|108|l . we 
find 



Q 



[n + l)"a' 



n „,2n 



From Eqs. 1114() and H115|l . we get 

(n-h l)"/3i/2' 



n z 



l/ri 



e{a) 



PGM 
R ' 



(115) 



(116) 



Eliminating M thanks to Eq. Hlll|) and introducing the 
Milne variables (|105|l . we finally obtain 



X 



Gi?2 



(n-f^ l)"/3i/2 



{uovn) — {l-vo).{in) 



For n — > -|-(xi, this parameter becomes equivalent to the 
one introduced in Sec. 12.41 for isothermal spheres (note 
that for sufficiently small index n, the parameter x can 
become negative when wq > 1; therefore, the definition of 
z given previously is not valid for small n). According to 
the turning point criterion, the series of equilibria becomes 
unstable in the grand canonical ensemble when dx/da = 
0. This gives the condition 



2(1 - Wo) - (n - l)uoWo = 0. 



(118) 



The values of a at which x is extremum can be determined 
by a graphical construction like in Chavanis (2002b). They 
correspond to the intersections between the hyperbole de- 
fined by Eq. (|118|l and the solution curve in the (u, v) 
plane (which depends on the index n). From these graph- 
ical constructions, we deduce that the function x{(^) h^-s 
one maximum for n < 5. For n = 5, the Lane-Emden 
equation can be solved analytically and we get 



X 



(l-l-ia2) 



12 ■ 



(119) 



which is maximum for a — (3/5)^^2, Por n > 5, the func- 
tion x(q^) an infinite number of extrema. The onset 
of gravitational instability corresponds to the first max- 
imum. The fugacity vs mass plot (for a given tempera- 
ture) can be deduced from Fig. El Collecting the results 
obtained by different authors, incomplete polytropcs with 



high density contrasts are unstable in MCE for n > 5, in 
CE for 71 > 3 and in GCE for all n. Complete polytropes 
are always stable in MCE, stable in CE only for 71 < 3 
and never stable in GCE. 

In the following, we analyze the stability of polytropic 
spheres by explicitly calculating the second order varia- 
tions of the thermodynamical potential. This study has 
been done by Taruya & Sakagami (2002a, b) in MCE and 
CE. We shall consider CE in relation with our previous 
study (Chavanis 2002b), GCE and GMCE. 

5.5. Generalized thermodynamical stability in the 
canonical ensemble 

In the canonical ensemble, the thermodynamical param- 
eter is the free energy (Massieu function) J = S — -^E 
which is explicitly given by 



J = I pd^r —I p^d^r, 



T 



2T 



(120) 



where we have used Eqs. (llOOfl and (|101|l . Using the equa- 
tion of state (|^ . the second order variations of free energy 
can be put in the form 

s'J = -^{J ^;^(^'^)'^'^ + / sps^d^y (121) 

As in previous papers, we introduce the function g(r) de- 
fined by 



Sp^ 



1 dq 



47rr2 dr ' 

and satisfying the perturbed Gauss theorem 

dS^ _ Gq 
dr r2 



(122) 



(123) 



The conservation of mass imposes q{0) = q{R) = 0. After 
an integration by parts, we can rewrite Eq. H121|) in the 
form 



2T 



drq 



dr \ Anr'^ dr 



G_ 



(124) 



The point of marginal stability is therefore determined by 
the condition 



d /p'^-2 



dr \ Aur"^ dr 



0. 



(125) 



This is the same equation as that determining the dy- 
namical stability of gaseous polytropes (Chavanis 2002b). 
We have solved this equation and showed that the bound- 
ary condition q{R) = returns the condition of Eq. (|112|l 
brought by the turning point criterion. The structure of 
the density perturbation profile that triggers instability is 
also discussed in our previous paper. 
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5.6. Generalized thermodynamical stability in the 
grand canonical ensemble 

The thermodynamical potential in the grand canonical en- 
semble is G = S" - + yM, i.e. 



G = nz'^/'' I pd^r - ^ / pd^r - ^ / p<i>dh. (126) 



Its second order variations are given by Eq. 11211) like for 
the free energy. Only the boundary conditions are differ- 
ent since the mass is not conserved anymore. Introducing 
the function X{r) defined by Eq. (|31|l . we obtain after an 
integration by parts 



6'G 



drX' 



G + Kj 



dr V 47r dr 



X'. (127) 



The point of marginal stability is determined by the equa- 
tion 



^ dr\ 47r 



GX' = 0, 



(128) 



with the boundary conditions (|32|l . Integrating once, we 
get 



K-i^ X" 

' At: 



GX = 0. 



In terms of the dimensionless variables introduced in Sec. 
15.31 we have to solve the problem 



1 



X = 0, 



X{0) = X'{a) = 0. 

Let us introduce the differential operator 



1. 



(130) 
(131) 

(132) 



1 

n a4 

Using the Lane-Emden equation H103|l . it is easy to check 
that 



>c(0 = e, C{C 



cm 



1 



£,9. (133) 



Therefore, the general solution of Eq. (|130|l satisfying 
X{0) = is 



n — \ 



(134) 



The critical value of a at which the series of equilibria 
becomes unstable in GCE is determined by the boundary 
condition X'{a) = 0. Using Eq. (|134|l and introducing the 
Milne variables (|l()5|l . we find that this condition is equiv- 
alent to Eq. H118|) in agreement with the turning point 
criterion. The density perturbation profile triggering in- 
stability can be deduced from Eqs. (jlTjl and ()134|l . We 
get 



5p 



in 



ci{vs - v), 



with Vs — 2/{n — 1). Equation (|135|l coincides with the 
expression found in CE. The number of nodes of the per- 
turbation profile can be determined by a graphical con- 
struction like in Chavanis (2002b). First, consider the case 
rt = 5 which can be solved analytically. The solution of the 
equation v{£,) = Vs is ^(i) = \/3. Since ^(i) > a = (3/5)^/^, 
we conclude that the perturbation profile in GCE has no 
node contrary to the equivalent situation in CE. This con- 
clusion remains true for all index n < 5. Indeed, the crit- 
ical case ^(1) = a would occur for an index n such that 
{us,Vs) belongs to the solution curve (substitute vq = Vs 
in Eq. (|118l) '). This is never the case for n < 5 (see 
Chandrasekhar 1942). Therefore, the structure of the den- 
sity profile for n < 5 is the same as for n = 5 so it has no 
node. For n > 5, the first mode of instability has no node, 
the second mode of instability one node etc... 

5. 7. Generalized thermodynamical stability in the 
grand microcanonical ensemble 

In the grand microcanonical ensemble, the fixed parame- 
ters are the fugacity and the energy. The normalized en- 
ergy is given by (Taruya & Sakagami 2002a) 



(129) A 



ER 

' GAP 



1 



n — 5 



1 



1 

vo 



2 Wo 



1 Vo 



(136) 



Eliminating f3 in Eq. H117() in profit of E, using Eqs. (|lll|l 
and H136|l . we find that the relevant control parameter 



{n 2l/«)2(2n-3) 



1 



{n + 1)4" Q4( 
is related to a by 



— —tG^rJe. 

n — 1 1 ' 



(137) 



71 — 2 Mo 

n + 1 uo 



(138) 



In the grand microcanonical ensemble, we have to max- 
imize the thermodynamical potential K. = S+^M at fixed 
E and p/T. The entropy and the energy are expressed in 
terms of p and K by Eqs. (|100|l . (|101|l and We vary 
the density and use the energy constraint to determine 
the corresponding variation of K (which is related to the 
inverse temperature /3). After some algebra, we find that 
the second order variations of JC are given by 



lJsp6<i>dh^\^J^{Sprd'r 



2n 



3(2n- 3) / p dH 



1 2 



(139) 



(To simplify the formulae, we have not written the positive 
proportionality factor in front of 6^IC). Introducing the 
variable X defined previously and integrating by parts, 
we can put the problem in the form 



(135) 6^)C = 



R 

dr / dr'X'{r)K{r,r')X'{r'), 
Jo 



(140) 
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with 



K{r,r') = |g + 7 



d f p d 
dr \Attp'^ dr 



J p d-^r 



3jp\ 
2p J _ 



(r) 



37P 



An 



2p 



3(2n- 3) 
(r'). (141) 



The condition of marginal stability reads 



GX 



IP 



X" 



4n 



1 



3(2n -3) Jp d^r 



37p\ 



2p 



(142) 



We now need to relate the gravitational potential $ to 
the Lane-Emden function 9. Integrating the equation of 
hydrostatic equilibrium for a polytropic equation of state, 
we readily find that 



1)^ = -$- 
P 



Const. 



(143) 



The constant can be determined from the boundary con- 
dition $(i?) = -GM/R. Using the relation 



GM 



-{n + l)Kpy''a6'{a), 



(144) 



which results from Eqs. (|110|) and H1U4|) . we find that Eq. 
p43fl can be written 



$ = (n + l)Kpl''\9{a) -9 + a9'{a)). 



(145) 



Inserting this relation in Eq. (|142|) and introducing the 
dimensionless variables defined in Sec. 15.31 we find that 



X 



_0l-nx" = 

n 

2(n + l)e 
J,"9^+H^d^ 
1 

^3 



(«) + 



3-2n. 



^9''XdC - 



2n 

a9'{a) + 9{a) 
2n-3 



a9\a) 
X{a) 



(146) 



where we have integrated by parts and used the iden- 
tity {^9)" = —^6'" equivalent to the Lane-Emden equation 
(|103() . We are led therefore to solve a problem of the form 



1 



X ^V9^ + W^, 



X{0) = X'{a) = 0. 



(147) 
(148) 



Using the identities l|133(l , we find that the general solution 
of Eq. lfTT7|l satisfying X{0) = is 



X = W^ + b 



^9+-{n-l)e0' 



2ai 



-Vi^9 



(149) 



The constant h is determined by the boundary condition 
X'{a) ~ 0. Substituting the resulting expression for X{£) 
in Eq. H146I) . and using the identities of Appendix IbI we 
obtain after tedious algebra a condition for a which is 
equivalent to the one obtained from the turning point ar- 
gument, i.e. by setting dv/da = 0. 



6. The isobaric ensemble 

We now consider the stability of isothermal spheres in 
contact with a medium exerting a constant pressure P 
on their boundary. In this isobaric ensemble, the volume 
of the system can change unlike in the Antonov prob- 
lem. This situation was first considered by Ebert (1955), 
Bonnor (1956) and McCrea (1957) who found that an 
isothermal gas becomes gravitationally unstable above a 
critical pressure Pmax or above a maximum compression. 
We shall complete these results by supplying an entirely 
analytical derivation of the stability criterion and by show- 
ing the equivalence between thermodynamical and dy- 
namical stability (Jeans problem). We shall also consider 
the case of polytropic gas spheres under external pressure. 

6.1. Isothermal gas spheres under external pressure 

Let us consider an isothermal gas sphere with mass M at 
temperature T and under pressure P. The first principle 
of thermodynamics can be written SS — ^SE — ^5V = 0. 
In the isothermal- isobaric ensemble (TEE), the relevant 
thermodynamical potential is the Gibbs energy Q — S — 
y-E—-^. Its first variations satisfy (5^ = —E5[^)—V5{y) 
so that, at statistical equilibrium, the system is in the state 
that maximizes Q at fixed temperature T and pressure P. 

For an isothermal gas, p = ^pT. Using Eq. ((T^ . the 
pressure on the boundary of the sphere is 



m p 



(150) 



The central density po can be eliminated via the relation 
(|15|l . Eliminating the radius R via the relation (|2()(l and 
introducing the Milne variables H16|l . we obtain after sim- 
plification 



V = AirG^P^A'PP 



UqVq. 



(151) 



Similarly, the normalized volume of the sphere is given by 
3V 1 



V 



47r/33G'3M3 ' 
Using the expansion of the Milne variables for a — > 0, 



(152) 



Q ^2 

uq = 6 a 

5 



(153) 



we find that the first order correction to Boyle's law due 
to self-gravity is 



PV = NkT - A 



GAP 
17W' 



(154) 



with A — i(|7r)^/'^. This formula can also be derived from 
the Virial theorem 



PV = I p d^r+l-W, 



(155) 



by using the isothermal equation of state and by assuming 
that the density p is uniform (in which case the potential 
energy W = —2>GM'^ JbR). As reported by Bonnor (1956), 
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0.2 0.3 
3V/47cP'g'm' 

Fig. 16. The P — V diagram for isothermal spheres. 
The series of equilibria becomes unstable at TBE in the 
isothermal-isobaric ensemble. 



Eq. I|154|l was first suggested by Terletsky (1952) as a sim- 
ple attempt to take into account gravitational effects in 
the perfect gas law. It appears here as a systematic ex- 
pansion in terms of the concentration factor a. For large 
concentrations, Eq. H154(l is not valid anymore and Eqs. 
H151(I - (I152|I must be used instead. The P — V curve is 
plotted in Fig. El It presents a striking spiral behaviour 
towards the limit point {Vs,Ps) = (1/8,8) corresponding 
to the singular sphere. This is similar to the E — T curve 
and other curves studied in Sec. 12.41 This is interesting to 
note, for historical reasons, because this diagram was dis- 
coved before the fundamental papers of Antonov (1962) 
and Lynden-Bell & Wood (1968) on the gravothermal 
catastrophe. Like the minimum energy of Antonov, there 
exists a maximum pressure Pmax = 1-40 k^T'^ / G'^ A'P 
above which no hydrostatic equilibrium for an isothermal 
gas is possible. Alternatively, this corresponds to a min- 
imum temperature Tmin = 0.919 G^/'^mM'^/'^P'^/'^/k for 
a given mass and pressure. From the turning point crite- 
rion, the series of equilibria becomes unstable in TBE at 
the point of maximum pressure. The condition dP/da = 
leads to 



(156) 

intersection 
the straight 
correspond- 
compressed 

0.29 G^M^rrv' /k^T^ are thermodynamically unstable. 

Using Eqs. ltTFT|) - (fTK^ and the identities l(T7J-|(lHl), we 
readily obtain 



2uo — uo = 0. 

The critical value a is determined by the 
between the spiral in the (u, v) plane and 
line H156|) . Its numerical value is a = 6.45 
ing to a density contrast of 14.0. More 
isothermal spheres, ie. those for which V 



dlnV 

da 



-("0 - 1), 



(157) 



and 

d\nP 1 



da 



(2uo - vq). 



Therefore, the "compressibility" k = 
can be expressed as 

1 2uo - vo 

K — . 

3 uo — I 



(158) 

-{d\nP/d\nV)N,T 
(159) 



Using uo = ZPV/NkT and wq = {At: /if/^GMm/kTV^/^ 
according to Eqs. I|151|l and (|152|l . we can check that Eq. 
((T^ is equivalent to Eq. (2.16) of Bonnor (1956), but it 
has been obtained here in a much more direct manner. 
Gravitational instability in TBE occurs when the com- 
pressibility becomes negative passing by k = 0. 

We now consider the dynamical stability of isothermal 
gas spheres under external pressure and show the equiv- 
alence with the thermodynamical criterion. This problem 
was studied numerically by Yabushita (1968) but we pro- 
vide here an entirely analytical solution. In a preceding 
paper, we have found that the equation of radial pulsa- 
tions for a gas with an isothermal equation of state could 
be written (Chavanis 2002a) 



1 dq 
m dr \ Ai:pr'^ dr 



Gq 



A2 



AnpTr 



(160) 



where A is the growth rate of the perturbation (such that 
dp ^ e'*'*) and q is defined by Eq. H122|) . In the isobaric 
ensemble, the boundary conditions are given by Eq. (jC.4|l 
of Appendix [CI Introducing dimensionless variables and 
considering the case of marginal stability A = 0, we are 
led to solve the problem 



d /e"^ dF 

d^[e^ 



= 0, 



with F{0) = and 
dF 



d^ d^ 



F = 0, 



at 



(161) 



(162) 



Now, the function F{£^) can be expressed as (Chavanis 
2002a) 



F{o^c,{ee-^~e^')- 



(163) 



Substituting this relation in the boundary condition 1162|) , 
using the Emden equation (|14l) and introducing the Milne 
variables (|16|) . we obtain the condition 



2ur) 



Vo = 0, 



(164) 



which coincides with the criterion (|156|l of thermodynam- 
ical stability. The expressions of the density perturbation 
profile and of the velocity profile at the point of marginal 
stability are given in Chavanis (2002a). Only the value of 
the critical point a changes. Due to the different bound- 
ary conditions, we now have Sv at the surface of the 
gaseous sphere since the sphere can oscillate. The per- 
turbation profile Sp/ p has one node for the first mode of 
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Fig. 17. Temperature vs enthalpy plot at fixed pressure 
and mass. The series of equilibria becomes unstable at 
HBE in the isenthalpic-isobaric ensemble. 



instability, two nodes for the second mode etc... The ve- 
locity profile Sv has no node for the first mode (except at 
the origin), one node for the second mode etc... 

We could also consider the situation in which the pres- 
sure P and the enthalpy H = E + PV are fixed (instead 
of the temperature). This corresponds to the isenthalpic- 
isobaric ensemble (HBE). Since 5S — ^5H — ^5P, the 
entropy must be a maximum at fixed H and P at statisti- 
cal equilibrium. Eliminating the radius R in Eqs. H19|l and 
in profit of the pressure P, we get after some algebra 



A' = 



H 



1 



(47r)l/4G3/4j(^3/2pl/4 (uoi;3)l/4 



7W0 



(165) 
(166) 



The T — H curve in HBE is plotted in Fig. [T7| It spirals 
around the fixed point (A'^,77;,) = (2i/'*/12, 23/-*) corre- 
sponding to the singular solution. There is no equilibrium 
state for A' > 0.262. This implies a minimum enthalpy 
for a fixed pressure or a maximum pressure for a fixed 
enthalpy. This is similar to the Antonov (1962) criterion, 
the pressure playing the same role as the box radius. The 
series of equilibria becomes unstable in HBE at the point 
of minimum enthalpy. The condition dK' /da ~ gives 



16uo + Swof — 38uo + 3wo = 0. 



(167) 



The critical value of a is 25.8 corresponding to a density 
contrast of 390. 

The specific heat at fixed pressure Cp ~ {dH / dT)N,p 
can be deduced from Eqs. p65|l and H166(l . After simplifi- 
cation, we find 



2 



Suqw ~ 38uo + Sfo 
2uo - Wo 



(168) 



For the ideal gas without gravity {uq,vo) = (3,0) we re- 
cover the well-known result Cp = ^Nk. Gravitational in- 
stability occurs in the isenthalpic-isobaric ensemble when 
Cp = passing from negative to positive values and in 
the isothermal-isobaric ensemble when Cp < 0. 



6.2. Polytropic gas spheres under external pressure 

We now consider a polytropic gas sphere of index n (or 
7 = 1 -|- 1/n) and mass M under pressure P. Using Eqs. 

and (|l()2|l , the pressure on the boundary of the sphere 
is given by 



P = Kp(Ry' = Kpl'^-eia 



(169) 



The central density po can be eliminated via the relation 
(|104|l . Eliminating the radius R via the relation (|110|l and 
introducing the Milne variables H105|l . we obtain after sim- 
plification 



K 



^-3 



P 



ri + l 
ri-3 



(170) 



A'3(n + l)-\ 

Similarly, the normalized volume of the sphere is given by 



V = W 



AttG^M'- 



1 



(woVq)"-' 



Using the expansion of the Milne variables for a 

1 o 



wn = 3 a 

5 



(171) 
0, 

(172) 



we find that the first order correction to the classical adi- 
abatic law due to self-gravity is 

CAP 



PV^ = KA'P - X- 



V— 



(173) 



Equation (|173|l can also be obtained from the Virial theo- 
rem (|155|) as before. The P — V curve is plotted in Figs. 1181 
and 1191 for different values of the polytropic index n. The 
volume is extremum at points a such that dV/da = 0. 
This leads to the condition uq — Ug. The pressure is ex- 
tremum at points a such that dV /da = 0. This leads to 
the condition 



2uq + (3 — n)uo = 0. 



(174) 



The number of extrema can be obtained by a graphical 
construction in the {u,v) plane like in Chavanis (2002b). 
For n < 3, there is no extremum for V and V. The pressure 
is a decreasing function of the volume and vanishes for 
a maximum volume Vmax corresponding to the complete 
polytrope (such that a — ^1). For n — 1, the solution of 
the Lane-Emden equation is = sin^/^ for £, < £,1 — tt 
and we explicitly obtain 



V 



(sin a — a cos a)^ 



V 



(175) 



for a < n. For 3 < n < 5, there is one extremum for V and 
V. This implies the existence of a maximum pressure Pmax 
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Using Eqs. l(T7n|l . l(T7T|) and the identities ((TnS|l and 
(I107|l . we readily obtain 



dlnV 

da 

and 

rflnT' 

da 



3 /n- 1 



a\n — 3 



Mo - 1 



{n — 3)a 



((3 - n)vQ + 2ua). 



(177) 



(178) 



Therefore, the "compressibility" k — {d\nP/d\nV)N,K 
can be expressed as 



1 / , 1 X 2uo + (3 - n)vo 
3 (n — l)uo — n + 3 



(179) 



10 



V We now consider the dynamical stability of polytropic 

gas spheres under external pressure. The equation of radial 
Fig. 18. The P-V diagram for polytropes with 1 < n < 3 pulsations for a gas with a polytropic equation of state can 



(specifically n — I) and 3 < n < 5 (specifically n = 4). 
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Fig. 19. The P — V diagram for polytropes with n > 5 
(specifically n — 10). 



and a minimum volume Vmm- In particular, for n = 5, the 
Lane-Emden equation can be solved analytically and we 
get 



V = 



36(1 



ia2)i2' 



V = 



a 



3^ 

15 



1 . 
1 + -.^ 



(176) 



The pressure is maximum for a — 3 and the volume is 
minimum for a — vTs. For n > 5 there is an infinity of 
extrema for V and V. The P~V curve spirals towards the 
limit point {Vs,Vs) corresponding to the singular sphere. 
From turning point arguments, the series of equilibria be- 
comes thermodynamically unstable at the point of max- 
imum pressure (for n > 3). Polytropes with n < 3 are 
always stable in the isobaric ensemble. 



be written (Chavanis 2002b) 



d f p-^-^ dq 
dr \ Anr-^ dr 



Gq 



A2 



(180) 



In the isobaric ensemble, the boundary conditions are 
given by Eq. HC.4|I of Appendix^ Introducing dimension- 
less variables and considering the case of marginal stability 
A = 0, we are led to solve the problem 



d_ 



g 1 — n 



dF 



nFiO 



with F(0) = and 
,dF d9 ^ „ 



0, 



(181) 



(182) 



Now, the function F(^) can be expressed as (Chavanis 
2002b) 



71 — 3 



2/)/ 



(183) 



Substituting this relation in the boundary condition H182|) , 
using the Lane-Emden equation H1U3|) and introducing the 
Milne variables H1(J5|) . we obtain the condition 



2uo + (3 — 71)110 = 0, 



(184) 



which coincides with the criterion (|174|l of thermodynam- 
ical stability. The expressions of the density perturbation 
profile and of the velocity profile at the point of marginal 
stability are given in Chavanis (2002b). Only the value of 
the critical point a changes. For 3 < 77 < 5, the density 
perturbation profile has one node and the velocity profile 
no node (except at the origin). For 77 > 5, the perturbation 
profile Sp/ p has one node for the first mode of instability, 
two nodes for the second mode etc... The velocity profile 
dv has no node for the first mode (except at the origin), 
one node for the second mode etc... 
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7. Stability of galaxies and stars 

In this section, we discuss in detail the relation between 
the dynamical and the thermodynamical stability of stel- 
lar systems and gaseous stars following our previous inves- 
tigations (Chavanis 2002a, b). This discussion should clar- 
ify the notion of generalized thermodynamics introduced 
by Tsallis (1988). 



7.1. Collisional vs collisionless relaxation 

Basically, a stellar system is a collection of N stars in 
gravitational interaction. Due to the development of en- 
counters between stars, this system is expected to relax to- 
wards a statistical equilibrium state. The relaxation time 
can be estimated by treiax ~ {N/ lnN)t]j where to is 
the dynamical time, i.e. the time it takes a star to cross 
the system (Binney & Tremaine 1987). For globular clus- 
ters, N ~ 10^, to ^ 10^ years so that the relaxation 
time is of the same order as their age (^ 10^° years). 
Therefore, globular clusters form a discret Hamiltonian 
system which can be studied by ordinary methods of sta- 
tistical mechanics. In particular, it can be shown rigor- 
ously that the Boltzmann entropy S'b[/] is the correct 
form of entropy for these systems (in a suitable ther- 
modynamic limit) although they are non-extensive and 
non-additive. The absence of entropy maximum in an un- 
bounded domain simply reflects the natural tendency of a 
stellar system to evaporate under the effect of encounters. 
However, evaporation is a slow process so that a glob- 
ular cluster passes by a succession of quasi-equilibrium 
states corresponding to truncated isothermals (see Sec. 
01). It should be emphasized that statistical mechanics is 
based on the postulate that the evolution is ergodic. Now, 
the kinetic theory of self-gravitating systems is not well- 
understood (Kandrup 1981). Indeed, it is possible to prove 
a H-theorem for the Boltzmann entropy only if simplify- 
ing assumptions are implemented (Markovian approxima- 
tion, local approximation, regularization of logarithmic di- 
vergences,...). It would be of interest to determine more 
precisely the influence of non-ideal effects on the relax- 
ation of stars. However, these effects are not expected to 
induce a strong deviation to Boltzmann's law. In particu- 
lar, the diffusion of stars is normal or slightly anomalous, 
i.e. polluted by logarithmic corrections (Lee 1968). Hence, 
Tsallis entropies, which often arise in the case of anoma- 
lous diffusion, do not seem to be justified in the context of 
"collisional" stellar systems. However, there are many ef- 
fects that can induce a lack of ergodicity and a deviation 
to Boltzmann's law. To a large extent, the dynamics of 
self-gravitating systems remains misunderstood and addi- 
tional work (using numerical simulations) is necessary to 
clarify the process of collisional relaxation. 

For elliptical galaxies, on the other hand, the relax- 
ation time due to close encounters is larger than their age 
by about ten orders of magnitude {N ^ 10^^, <_d ^ 10^ 
years, age ~ 10^° years). Therefore, their dynamics is es- 
sentially encounterless and described by the self-consistent 



Vlasov-Poisson system. Elliptical galaxies form therefore 
a continuous Hamiltonian system. Now, due to a compli- 
cated mixing process in phase space, the Vlasov-Poisson 
system can achieve a metaequilibrium state (on a coarse- 
grained scale) on a very short timescale. A statistical the- 
ory appropriate to this violent relaxation was developed 
by Lynden-Bell in 1967 (see a short review in Chavanis 
2002f). For collisionless stellar systems, the infinite mass 
problem is solved by realizing that violent relaxation is in- 
complete so that the galaxies are more confined that one 
would expect on the basis of statistical considerations. 

Collisionless and collisional relaxations are therefore 
two distinct processes in the life of a self-gravitating sys- 
tem and they must be studied by different types of statis- 
tical mechanics. Note that a globular cluster first experi- 
ences a collisionless relaxation leading to a virialized state 
and then evolves more slowly towards another equilibrium 
state due to stellar encounters. We can discuss these two 
successive equilibrium states in a more formal manner. Let 
us consider a fixed interval of time t and let the number of 
particles N +oo. In that case, the system is rigorously 
described by the Vlasov equation, and a metaequilibrium 
state is achieved on a timescale independant on iV as a 
result of violent relaxation. Alternatively, if we fix N and 
let t +00, the system is expected to relax to an equi- 
librium state resulting from a collisional evolution. These 
two equilibrium states are of course physically distinct. 
This implies that the order of the limits N — > -|-oo and 
t +00 is not interchangeable for self-gravitating sys- 
tems (Chavanis 2002e). 

7.2. The Vlasov-Poisson system 

In Sees. 12141 we have implicitly considered the case of col- 
lisional stellar systems, such as globular clusters, relax- 
ing via two-body encounters. We now turn to the case 
of collisionless stellar systems, such as elliptical galaxies, 
described by the Vlasov equation 



df df df 

— -f V— + F— = 0, 

dt dr dv ' 



(185) 



coupled with the Poisson equation |(2Jl. The Vlasov equa- 
tion (|185|l simply states that, in the absence of encoun- 
ters, the distribution function / is conserved by the flow 
in phase space. This equation conserves the usual invari- 
ants (0) and (O but also an additional class of invariants 
called the Casimirs 



h = J h{f)d^rd^v, 



(186) 



where h is any continuous function of /. This infinite set 
of additional constraints is due to the collisionless nature 
of the evolution and results from the conservation of / 
and the incompressibility of the flow in phase space. As 
we shall see, they play a crucial role in the statistical me- 
chanics of violent relaxation. We note that the conserva- 
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tion of the Casimirs is equivalent to the conservation of 
the moments 

Mn = J /"d^rd^v. (187) 
It is well-known that any distribution function of the 

2 

form / — /(e), where e = ^ -f is the individual energy 
of a star, is a stationary solution of the Vlasov equation. 
This is a particular case of the Jeans theorem (Binney 
& Tremaine 1987). We shall assume furthermore that 
/'(e) < so that high energy stars are less probable than 
low energy stars. Any such distribution function can be 
obtained by extremizing a functional of the form 

S[f] = -Jc{f)d^rd^v, (188) 

at fixed mass M and energy E, where C(/) is a convex 
function, i.e. C"{f) > 0. Furthermore, maxima of S[f] at 
fixed mass and energy are dynamically stable with respect 
to the Vlasov-Poisson system (Ipser 1974, Ipser & Horwitz 
1979). This is a sufficient condition of nonlinear stabil- 
ity. We are tempted to believe that it is also a necessary 
condition of nonlinear stability although this seems to go 
against the theorem of Doremus et al (1971) which states 
that all distribution functions with /'(e) < are dynam- 
ically stable (even if they are minima or saddle points of 
S at fixed E and M). In fact, this criterion is only a con- 
dition of linear stability (Binney & Tremaine 1987) and 
it may not be valid for box confined models (its general 
applicability has also been criticized). In the following, we 
shall consider that / = /(e) is dynamically stable "in a 
strong sense" if and only if it is a maximum of S[f] at fixed 
E and M. This is equivalent to / = /(e) being a minimum 
of E at fixed S and M (Ipser & Horwitz 1979). It can be 
shown that stellar systems satisfying this requirement are 
necessarily spherically symmetric. 

Introducing appropriate Lagrange multipliers and 
writing the variational principle in the form 

SS - (55E - a5M = 0, (189) 

one finds that 

C'(/) = -/3e - a. (190) 

Since C" is a monotonically increasing function of /, we 
can inverse this relation to obtain 

J^F{Pe + a), (191) 

where F{x) = {C')^^{—x). From the identity 

/'(e) = -/3/C"(/), (192) 

resulting from Eq. p9()(l . /(e) is a monotonically decreas- 
ing function of energy if /? > 0. 

We note also that for each stellar systems with / = 
/(e), there exist a corresponding barotropic gas with 
the same equilibrium density distribution. Indeed p = 



/ /d^v = p(<i>), P = I J fv'^d^v = p($), hence p = p{p). 
Writing explicitly the pressure in the form 

P=l f F(/3e-ha)47r[2(e-$)]3/2de, (193) 

and taking its gradient, we obtain the condition of hydro- 
static equilibrium 

Vp = -pV$. (194) 

Therefore, the Jeans equations of stellar dynamics reduce 
in these cases of spherical systems with isotropic pressure 
to the equation of hydrostatics for a barotropic gas. 

7.3. Violent relaxation of collisionless stellar systems 

The Vlasov equation does not select a universal form of 
entropy S[f], unlike the Boltzmann equation for example. 
Therefore, the notion of equilibrium is relatively subtle. 
When coupled to the Poisson equation, the Vlasov equa- 
tion develops very complex filaments as a result of a mix- 
ing process in phase space. In this sense, the fine-grained 
distribution function /(r, v, t) will never converge towards 
a stationary solution /(r,v). This is consistent with the 
fact that the functionals S[f] are rigorously conserved by 
the fiow (they are particular Casimirs). However, if we 
introduce a coarse-graining procedure, the coarse-grained 
distribution function /(r, v, t) will relax towards a metae- 
quilibrium state /(r,v). This collisionless relaxation is 
called violent relaxation (Lynden-Bell 1967) or chaotic 
mixing. It can be shown that S[f] calculated with the 
coarse-grained distribution function at time f > is larger 
than S[f] calculated with the fine-grained distribution 
function at t = (a monotonic increase of H is not im- 
plied). This is similar to the i7-theorem in kinetic theory 
except that S[f] is not necessarily the Boltzmann entropy. 
For that reason, S[f] is sometimes called a H-function 
(Tremaine et al 1986). Boltzmann and Tsallis functionals 
Sb [/] and Sq [/] are particular H- functions corresponding 
to isothermal stellar systems and stellar polytropes. 

It is important to note, furthermore, that among the 
infinite set of integrals conserved by the Vlasov equation, 
some are more robust than others. For example, the mo- 
ments AI^^-^-[f] ~ J f d'^rd'^v calculated with the coarse- 
grained distribution function vary with time for n > 1 
since f ■ By contrast, the mass and the energy 

are approximately conserved on the coarse-grained scale. 
We shall therefore classify the collisionless invariants in 
two groups: the energy E and the mass M will be called 
robust integrals since they are conserved by the coarse- 
grained dynamics while the moments Af„ with n > 1 
will be called fragile integrals since they vary under the 
operation of coarse-graining. Note that, in practice, the 
moments are also altered at the fine-grained scale 

because, intrinsically, the system is not continuous but 
made of stars. Therefore, the graininess of the medium 
can break the conservation of some moments, in general 
those of high order, even if we are in a regime where the 
evolution is essentially encounter less. 
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7.4. Lynden-Bell's theory of violent relaxation 

Now, the question of fundamental interest is the following: 
given some initial condition, can we predict the equilib- 
rium distribution function that the system will eventually 
achieve? Lynden-Bell (1967) tried to answer this question 
by using statistical mechanics arguments. He introduced 
the density probability p(r, v, rj) of finding the value of the 
distribution function rj in (r, v) at equilibrium. Then, us- 
ing a combinatorial analysis he showed that the relevant 
measure is the Boltzmann entropy 



S[p] = - J phipd^rd^ 



drj, 



(195) 



which is a functional of p. This mixing entropy has been 
justified rigorously by Michel & Robert (1994), using the 
concept of Young measures. Therefore, although the sys- 
tem is non-extensive, the Boltzmann entropy is the correct 
measure. Assuming ergodicity (which may not be realized 
in practice, see Sec. I7.5|l the most probable equilibrium 
state, i.e. most mixed state, is obtained by maximizing 
S[p] taking into account all the constraints of the dynam- 
ics. The optimal equilibrium state can be written 



(196) 



p(r,v,,7) = |x(ry)e-('''+")^ 

where x(?7) = exp(— X)n>i o^nV") accounts for the conser- 
vation of the fragile moments M„>i = J prf-drfd^rd^-v and 
a, P are the usual Lagrange multipliers for M and E (ro- 
bust integrals). The distinction between these two types of 
constraints is important (see below) and was not explicitly 
made by Lynden-Bell (1967). The "partition function" 



r+oo 

Z = I xiv) e 
Jo 



dr], 



(197) 



is determined by the local normalization condition 
/ pdrj = 1 . The condition of thermodynamical stability 
reads 

(5V =-- [ '^^Pf-d^rd^wdri ~-(3 I S^Sjd^rd\ < 0, 
2 J p 2 J 

(198) 

for all variation Sp{r, v, rj) which does not change the con- 
straints to first order. The equilibrium coarse-grained dis- 
tribution function / = / pr]dri can be expressed as 

Taking the derivative of Eq. (|199|l , it is easy to show that 



(199) 



/ (e) = -/3/2 



/2 



p{-q~ ffdri > 0, 



(200) 



where /2 is the centered local variance of the distribu- 
tion p(r,v,77). The integrability condition J fd^v < oo 
requires that the temperature is positive {f3 > 0). Then, 
from Eq. 1200() . we deduce that / is a monotonically de- 
creasing function of e. Therefore, from Sec. 17.21 we know 
that for each function /(e) of the form p99|l there exists a 



H-f unction S[f] whose extremization at fixed mass and en- 
ergy returns Eq. (|199|l . Furthermore, it can be shown that 
the condition (|198|l of thermodynamical stability (max- 
imum of S[p] at fixed E and M„) implies that / is a 
nonlinearly stable solution of the Vlasov-Poisson system. 
Therefore, it is a maximum of S[f] at fixed mass and 
energy. This functional depends on the initial conditions 
through the Lagrange multipliers a„>i, or equivalently 
through the function xiv)! which are determined by the 
fragile moments Mn>i. Therefore, S[f] is non-universal 
and can take a wide variety of forms. In general, S[f] dif- 
fers from the ordinary Boltzmann entropy Ssif] derived 
in the context of coUisional stellar systems. This is due 
to the additional constraints (Casimir invariants) brought 
by the Vlasov equation. For a given initial condition, the 
H- function S[f] maximized by the system at statistical 
equilibrium is uniquely determined by the statistical the- 
ory according to Eqs. H190() and H199() . 



In general, S[f] is not an entropy, unlike S[p], as it 
cannot be directly obtained from a combinatorial analy- 
sis. There is a case, however, where S[f] is an entropy. 
If the initial condition is made of patches with value 
/o — Vo and fa — (two-levels approximation), then 
the entropy S[p] is equivalent to the Fermi-Dirac en- 
tropy Sp.D.if] = -/[//% In 7/r;o + (1 - 7/7?o)ln(l - 
f /riQ)]d^rd^v and Lynden-Bell's distribution 11991) co- 
incides with the Fermi-Dirac distribution (Lynden-Bell 
1967, Chavanis & Sommeria 1998). Moreover, all the mo- 
ments can be expressed in terms of the mass alone as 
Af„ = t]qM. In the non-degenerate (or dilute) limit / ^ 
770 , the Fermi-Dirac entropy reduces to the Boltzmann en- 
tropy 5's [/] = —/ / In fd^rd^v. Quite generally, in the di- 
lute limit, the equilibrium distribution is a sum of isother- 
mal distributions. 



We can note some general predictions of the statistical 
theory, which are valid for any initial condition /o(r, v). (i) 
The predicted distribution function is a function / = /(e) 
of the stellar energy alone in the non-rotating case and a 
function / = /(ej) of the Jacobi energy ej = e — • (r x v) 
alone in the rotating case, (ii) /'(e) < so that high en- 
ergy stars are less probable than low energy stars, (iii) 
/(e) is strictly positive for all energies (it never vanishes), 
(iv) /(e) is always smaller than the maximum value f^°^^ 
of the initial condition, (v) According to Eqs. H200|) and 
(|192|l . we have the general relation /2 — 1/C"{f). (vi) For 
a box-confined system, there always exist a global maxi- 
mum of entropy. Therefore, contrary to coUisional stellar 
systems, it is not necessary to introduce a small-scale cut- 
off to make the thermodynamical approach rigorous. This 
is due to the Liouvillc theorem which puts an upper bound 
on the value of the distribution function like the Pauli ex- 
clusion principle in quantum mechanics (see Chavanis et 
al 1996, Chavanis & Sommeria 1998). 
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7.5. Incomplete relaxation 

The prediction of Lynden-Bell's statistical theory crucially 
depends on the initial conditions. This is because we need 
to know the value of the moments M„>i which can only 
be deduced from the fine-grained distribution function at 
t = Q (say). Indeed, these moments are fragile, or mi- 
croscopic, constraints which cannot be mesured from the 
coarse-grained flow at t > since they are altered by the 
coarse-graining procedure. This is a specificity of contin- 
uous Hamiltonian systems which contrasts with ordinary 
(discrete) Hamiltonian systems for which the constraints 
are robust, or macroscopic, and can be evaluated at any 
time (as the mass M and the energy E). Now, Lynden- 
Bell (1967) gives arguments according to which elliptical 
galaxies should be non-degenerate (except possibly in the 
nucleus). In that case, the dilute limit of his theory applies 
in the main body of the galaxy and leads to an isothermal 
distribution. The theoretical justification of an isothermal 
distribution for coUisionless stellar systems was considered 
as a triumph in the 1960's. Indeed, Lynden-Bell's statisti- 
cal theory of violent relaxation could explain the observed 
isothermal core of elliptical galaxies without recourse to 
collisions that operate on a much longer timescale. 

However, it is easy to see that the statistical theory 
of violent relaxation cannot describe the halo of elliptical 
galaxies. Indeed, at large distances, the density decreases 
as resulting in the infinite mass problem. Lynden- 
Bell (1967) has related this mathematical difficulty to the 
physical problem of incomplete relaxation and argued that 
his distribution function has to be modified at high en- 
ergies. Indeed, the relaxation is effective only in a finite 
region of space and persists only for a finite period of 
time. Therefore, there is no reason to maximize entropy 
in the whole available space. The ergodic hypothesis which 
sustains the statistical theory applies only in a restricted 
domain of space, in a sort of "maximum entropy bubble" , 
surrounded by an unmixed region which is only poorly 
sampled by the system. Accordingly, the physical picture 
that emerges is the following: during violent relaxation, 
the system has the tendency to reach the most mixed 
state described by the distribution 1)199(1 . However, as it 
approaches equilibrium, the mixing becomes less and less 
efficient and the system settles on a state which is not the 
most mixed state. 

How can we take into account incomplete relaxation? 
A first possibility is to confine the system artificially in- 
side a box, where the box delimits the typical region in 
which the statistical theory applies (Chavanis & Sommeria 
1998). Another possibility is to use a parametrization of 
the coarse-grained dynamics of coUisionless stellar systems 
in the form of relaxation equations involving a space and 
time dependant diffusion coefficient related to the fine- 
grained fluctuations of the distribution function (Chavanis 
et al 1996, Chavanis 1998). This space and time dependant 
diffusion coefficient can freeze the system in a "maximum 
entropy bubble" and lead to the formation of self-confined 
stellar systems. Finally, a third possibility is to construct 



a model of incomplete violent relaxation which provides a 
depletion of stars in the high energy tail of the distribu- 
tion function. If the system is isolated, as most elliptical 
galaxies are, Hjorth & Madsen (1993) have proposed a 
dynamical scenario attempting to take into account in- 
complete relaxation in the halo. They introduce a two- 
step process: (i) in a first step, they assume that violent 
relaxation proceeds to completion in a finite spatial re- 
gion, of radius rmax, which represents roughly the core of 
the galaxy (where the fluctuations are important). At this 
stage, the escape energy represents no special threshold 
so that negative as well as positive energy states are pop- 
ulated in that region, (ii) After the relaxation process is 
over, positive energy particles leave the system and par- 
ticles with '^{rmax) < e < move in orbits beyond r^nax, 
thereby changing the distribution function to something 
significantly 'thinner' than a Boltzmann distribution. The 
crucial point to realize is that the differential energy dis- 
tribution iV(e), where N{e)de is the number of stars with 
energy between e and e -I- de, will be discontinuous at 
the escape energy e = since there is a flnite number 
of particles within r„iax after the relaxation process. It 
can be shown that this discontinuity implies necessarily 
that /(e) - {-ef/"^ for e ^ 0" (Jaffe 1987). Therefore, 
the distribution function consistent with this scenario is 
(Hjorth & Madsen 1993) 

( Ae-f^"^"' $0 < e < e', 
f=\B (-e)5/2 e' < e < 0, (201) 
I e > 0. 

This model corresponds formally to a composite configu- 
ration with an isothermal core in which violent relaxation 
is efficient and a polytropic halo (with index n = 4) which 
is only partially relaxed (in reality, there is a small devia- 
tion to the ideal value n = 4 because the potential created 
by the core is not exactly Keplerian; see Hjorth & Madsen 
1993). For this model, the density decreases as r^^ like in 
elliptical galaxies. Furthermore, this model can reproduce 
de Vaucouleur's R^/'^ law for the surface brightness of el- 
lipticals (Binney & Tremaine 1987). Alternatively, if the 
system is subject to tidal forces (this may be the case for 
dark matter halos), one can propose an extended Michie- 
King model of the form (Chavanis 1998) 

'^o x+e-P^o^ £ < e™, (202) 

e > Em, 

which takes into account the specificities of Lynden-Bell's 
distribution function (in particular the exclusion principle 
associated with Liouville's theorem) and leads to a con- 
finement of the density at large distances (i.e. the density 
drops to zero at a finite distance identified with the tidal 
radius). 

7.6. Tsallis entropies 

The statistical theory of Lynden-Bell encounters some 
problems because the distribution function /(e) never van- 
ishes. Therefore, we have to invoke incomplete relaxation 
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and introduce dynamical constraints in order to explain 
the observation of confined galaxies and solve the infinite 
mass problem. This is probably the correct approach to 
the problem. Another line of thought, defended by Tsallis 
and co-workers, is to change the form of entropy in order 
to obtain a distribution function that vanishes at some 
maximum energy. Therefore, Tsallis generalized thermo- 
dynamics is essentially an attempt to take into account 
incomplete mixing and lack of ergodicity in complicated 
systems. If we assume that the system is classically de- 
scribed by the Boltzmann entropy Sslf] (which is true 
for coUisionless stellar systems only in a coarse-grained 
sense and in a dilute limit), then Tsallis proposes to re- 
place SbIJ] by a g-entropy Sq[f]. The parameter q mea- 
sures the efficiency of mixing (q = 1 ii the system mixes 
well) but is not determined by the theory. 

This approach is attractive but the justification given 
by Tsallis and co-workers is misleading. Fundamentally, 
ordinary thermodynamics (in the sense of Lynden-Bell ) 
is the correct approach to the problem. The Boltzmann 
entropy S[p] has a clear physical meaning as it is propor- 
tional to the logarithm of the disorder, where the disorder 
is equal to the number of microstates consistent with a 
given macrostate. Of course, if the system does not mix 
well, statistical mechanics loses its power of prediction. 
However, the state resulting from incomplete violent re- 
laxation is always a nonlinearly stable stationary solution 
of the Vlasov equation on a coarse-grained scale (i.e., for 
/). Therefore, if / = /(e), it maximizes a if-function 
S[f] = - J C{f)d^T(p-v at fixed mass M and energy E. 
The function C(/) is influenced by thermodynamics (mix- 
ing) but it cannot be determined by pure thermodynami- 
cal arguments because of incomplete relaxation. 

The iJ-function S[J] selected by the system depends 
both on the initial condition and on the type of mixing. 
Tsallis entropy Sq[f] is just a particular H- function lead- 
ing to confined structures (stellar polytropes). Now, it 
is well-known that stellar polytropes give a poor fit of 
elliptical galaxies. In particular, they cannot reproduce 
de Vaucouleur's law, neither the r""* density decrease of 
these objects. Indeed, for stellar polytropes, the density 
decreases as for n = 5, the density drops to zero at 
a finite distance for n < 5 and the mass is infinite for 
n > 5. We conclude, therefore, that observations of ellip- 
tical galaxies do not support the prediction of Tsallis. The 
model of Hjorth & Madsen (1993) based on the statistical 
theory of Lynden-Bell improved by a physically motivated 
scenario of incomplete relaxation provides a much better 
agreement with observations. This equilibrium state can- 
not be described by Tsallis distribution which fails to re- 
produce simultaneously an isothermal core and an exter- 
nal confinement. In fact, the model of Hjorth & Madsen 
(1993) can be obtained by maximizing a H-function with 

c(/)={^{;-{V)j;^:: « 



at fixed mass and energy. Therefore, their model is consis- 
tent with a composite Tsallis model with g = 1 (complete 
mixing) in the core and q = {2n — l)/(2n — 3) = 7/5 
(incomplete mixing) in the halo. More generally, we could 
argue that q varies in space in order to take into account a 
variable mixing efficiency (very similarly, stars are some- 
times approximated by composite polytropic models with 
a space dependent index in order to distinguish between 
convective and radiative regions). However, this gives to 
Tsallis distribution a practical, not a fundamental justi- 
fication. It can just fit locally the distribution function 
of coUisionless stellar systems resulting from incomplete 
violent relaxation. 

In addition, not all stable stationary solutions of the 
Vlasov equation can be obtained by maximizing a H- 
function at fixed mass and energy. This maximization 
problem leads to f = /(e) with /'(e) < 0. For spheri- 
cal stellar systems, the distribution function can depend 
on the specific angular momentum j = r x v in addi- 
tion to the stellar energy e and this is what happens 
in certain models of incomplete violent relaxation with 
anisotropic velocity distribution (Lynden-Bell 1967). On 
the other hand, for rotating stellar systems, the maximiza- 
tion of S[f] at fixed mass, energy and angular momentum 
predicts that / = /(ej) with /' < 0, where ej is the 
Jacobi energy. However, it is well-known that such distri- 
bution function cannot account for the triaxial structure 
of elliptical galaxies. More general distribution functions 
must be constructed in agreement with the Jeans theorem 
(Binney & Trcmaine 1987). However, it is possible that rel- 
evant distribution functions maximize a i?-function with 
additional constraints. These constraints may correspond 
to adiabatic invariants that are approximately conserved 
during violent relaxation. For example, anisotropic veloc- 
ity distributions can be obtained by maximizing a H- 
function while conserving the individual distributions of 
angular momentum L2 = / /(r x v)^d^rd'^v (or other dis- 
tribution) in addition to mass and energy. This leads to 
-2 

/ = /(ca) with ea = and /'(£„) < {va is called the 

anisotropy radius). It is expected that this maximization 
problem is a condition of nonlinear stability. For rotating 
systems, it is not known whether we can construct distri- 
bution functions yielding triaxial bodies by maximizing a 
if-function at fixed M, E, L and additional (adiabatic) 
invariants. These ideas may be a route to explore. 

7. 7. Dynamical stability of coUisionless stellar systems 

On a formal point of view, the maximization problem de- 
termining the nonlinear stability of a coUisionless stellar 
system is similar to the one which determines the thermo- 
dynamical stability of a coUisional steUar system. In ther- 
modynamics, we maximize the Boltzmann entropy <S'b[/] 
at fixed mass M and energy E in order to determine the 
most probable equilibrium state. In stellar dynamics, the 
maximization of a H-function S[f] at fixed mass and en- 
ergy determines a nonlinearly stable stationary solution of 
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the Vlasov-Poisson system. Therefore, we can use a ther- 
modynamical analogy to analyse the dynamical stability 
of collisionless stellar systems. In this analogy, the func- 
tional S[f] plays the role of a generalized entropy, the 
Lagrange multiplier /3 plays the role of a generalized in- 
verse temperature f3 = 1/T, the curve P{E) the role of a 
generalized caloric curve and the maximization of S[f] at 
fixed E and M corresponds to a microcanonicial descrip- 
tion. This analogy explains why we take C to be convex 
(and not concave) in Eq. H188(l . Then, we must look for 
maxima of S (not minima) like in thermodynamics. We 
can also introduce a J- function defined by J[/] = S* — f3E 
which is similar to a free energy in thermodynamics (J is 
the Legendre transform of S). The maximization of J[f] 
at fixed M and /3 corresponds to a canonical description. 

The condition that / is a maximum of S at fixed mass 
and energy is equivalent to the condition that S^J — S^S— 
PS'^E is negative for all perturbations that conserve mass 
and energy to first order. This condition can be written 



6' J 



1 



2 

ysf 



(3 / SpS<^d^r < 0, 



SE = SM = 0. 



(204) 



The condition that / is a maximum of J = S~ j3E at fixed 
mass and temperature is equivalent to the condition that 
5^, J is negative for all perturbations that conserve mass. 
This can be written 



5^ J 



C"{f)i^dhd''^^^l3 I SpS'Pdh < 0, 



V 6f I 6M = 0. 



According to Eq. (fTi^ . S'^J can be re-written 



6' J 



-J 



-d'^rd-^v 



SpS^d^r 



(205) 



(206) 



which is the usual form considered in the literature. We 
emphasize, however, that the role of the first order con- 
straints is crucial in the stability analysis and is often un- 
derestimated. If / is a maximum of J at fixed M and /3 
then it is necessarily a maximum of S at fixed M and E. 
Indeed, if the inequality 1)2051) is satisfied for all perturba- 
tions that conserve mass, it is a fortiori satisfied for pertur- 
bations that conserve mass and energy. However, the recip- 
rocal is wrong in general. In thermodynamics, this means 
that canonical stability implies microcanonical stability 
(but the converse is wrong in general). Exploiting this 
thermodynamical analogy, we can use the turning point 
arguments of Katz (1978) to settle whether a solution of 
the form (|190|l is a maximum or a minimum (or a saddle 
point) of S or J. The canonical and microcanonical de- 
scriptions will be inequivalent if the curve f3{E) presents 
turning points leading to regions where C — dE/dT < 0. 
This corresponds to regions of negative specific heats in 
thermodynamics. In these regions, the curve S{E) has a 
convex dip. This convex dip is the signal of phase transi- 
tions in thermodynamics. 



We can also try to reduce the stability problem to 
the study of an eigenvalue equation. To that purpose, it 
is convenient to first maximize S[f] at fixed mass, en- 
ergy and density p(r) in order to obtain a functional S[p]. 
Now, maximizing S[f] at fixed mass, energy and density 
is equivalent to maximizing S[f] at fixed kinetic energy K 
and density. Introducing a Lagrange multiplier /3 for the 
kinetic energy and a space dependant Lagrange multiplier 
A(r) for the density constraint, we get 

,,2 



f = F 



A(r) 



(207) 



where F is the same function as in Sec. 17.21 Since 
5'^S = - I C"(/)^^d^rd^v < 0, (208) 



we conclude that Eq. (|207() is a maximum of S at fixed 
K and p. This is also a maximum of J = S* — PE at fixed 
[3 and p. 

With the distribution function (|207|l . we can write the 
density and the pressure as 



1 



' /33/2 
where 



5(/3A), 



1 



r + oo 

g{x) ^ / 47rV2F(i + x)t^/^dt, 
Jo 



and 

h{x] 



-F{t + x)t^l^dt. 



(209) 



(210) 



(211) 



We can now express E and S as functionals of p and /3. 
The energy ^ is simply given by 



On the other hand, according to Eq. H188|) . we have 

//•+00 
d^r J C[F{Pe)]{e- Xy^^de. 



(212) 



(213) 



Setting t = (3{e — A), we get 
47rV2 '■+°° 



5 = 



^3/2 



d^'r 



C[F{t + I3X)] t^'^dt. 



(214) 



Integrating by parts and using C'[F{x)\ = — x, wc find 
that 



S = 



87rV2 
3/33/2 



F'{t + pX){t + f3X)t^/^dt. (215) 



Integrating by parts one more time and using Eqs. H209|) . 
(I210|l and (|211|l . we finally obtain 



S^-P I pd-'r + (3 I Xpd-'r. 



(216) 



We now consider the maximization of S[p] at fixed 
mass and energy. This is similar to a microcanonical situ- 
ation in thermodynamics. In that case, (3 is not fixed and 
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its variation 6/3 is related to the change of density Sp via 
the energy constraint (|212|l . using Eqs. H209() . piOl) and 
(|211|) . This implies in particular that the barotropic rela- 
tion p = p{p) , valid at equilibrium, is not preserved when 
we consider variations around equilibrium. Indeed, elimi- 
nating j3\ in Eq. we find that p = p-^/'^H{P^/'^p) 
so that (3 enters explicitly in the relation between p and p. 
In the case of isothermal stellar systems and stellar poly- 
tropes described by Boltzmann and Tsallis entropies, it 
is possible to reduce the stability problem to the study 
of an eigenvalue equation (Padmanabhan 1989, Taruya & 
Sakagami 2002a). However, it seem difficult to make this 
reduction in the general case. 

7.8. Dynamical stability of gaseous spheres 

We now consider the maximization of J[p] = S[p] — /3E[p] 
at fixed (3 and M . This is similar to a canonical situation 
in thermodynamics. According to Eqs. (|212|l and (|216|) . 
we have 



J = 13 pcPr + (3 



p^crr, 



(217) 



Since f3 is constant, the barotropic relation p = p(p) re- 
mains valid around equilibrium. Using the relation h'{x) = 
—g{x) obtained from Eqs. H210(l and H211|l by a simple in- 
tegration by parts, it is easy to check that Eq. (|209|l yields 



X + ^ = - 

P Jo P 



'pip'), , 
-^dp. 



(218) 



Therefore, the functional H217|) can be rewritten as 



J 



-P / P 



dp'd-'r - -(3 / p<^d^r 



(219) 



If we cancel the first order variations of J at fixed mass and 
f3, we get A = $ + a//3 and we recover Eq. H191|) . On the 
other hand, the condition that Eq. H191() is a maximum 
of J at fixed M and (3 requires that the second order 
variations of J be negative. This leads to the inequality 



-13 I ^-^{Spfd'^v - i/3 / SpS^d'^v < 0, (220) 



for all perturbations dp that conserve mass. Using the con- 
dition of hydrostatic equilibrium, (5^ J can be re- written 



1 



JSpf 

-p'i'^y 



5p6^d^r 



(221) 



We note, in passing, the similarity with the pseudo-energy 
(or Arnol'd invariant) in two-dimensional hydrodynam- 
ics (see Chavanis 2002h). Therefore, Arnol'd theorems 
can be directly extended to (3D) self-gravitating bod- 
ies. In particular, in a bounded domain, the system is 
nonlinearly stable, for any type of perturbation, if < 
—p'{^) ^ Ai/(47rG), where Ai is the smallest eigenvalue 
of —A (for isothermal spheres this corresponds to a < tt 
and for polytropic spheres to a < tt/^/ti where a is de- 
fined by Eqs. H15|l and 1104(1 respectively). Furthermore, 



we have already shown in the case of isothermal stellar 
systems and stellar polytropes that the stability condi- 
tion H220|l can be reduced to the study of an eigenvalue 
equation. This reduction can be generalized to any stellar 
system described by a H- function S[f] of the form H188|) . 
Considering spherically symmetrical perturbations, intro- 
ducing the variable q{r) defined by Eq. H122|l and repeat- 
ing exactly the same steps as for isothermal and polytropic 
spheres (see Chavanis 2002a and Sec. 15.511 . we end up on 
the eigenvalue problem 



d f P'jp) d 
dr \47rpr2 dr 



G 



qxir) = Xqx{r), 



(222) 



with appropritate boundary conditions (see Appendix lO). 
If all the A's are negative, then Eq. H191() is a maximum 
of J at fixed M and f3. 

We now consider the stability of a gaseous system with 
an equation of state p = p{p). The energy of this gaseous 
configuration is (Lynden-Bell & Sanitt 1969) 



W = 



(223) 



We observe that the J-function H219() associated with a 
stellar system is equal (up to a negative proportional- 
ity factor —(3) to the energy H223|l of the corresponding 
barotropic gas. The stability theorem of Chandrasekhar 
for gaseous spheres requires that S^W > (see Binney 
& Tremaine 1987). This leads to the inequality (|220|l 
and to the eigenvalue equation H222|l . We now consider 
the dynamical stability of a self-gravitating gaseous sys- 
tem with respect to the Navier-Stokes or Euler equations 
(Jeans problem). This problem can be solved explicitly 
for isothermal and polytropic spheres (Chavanis 2002a, b). 
Adapting the same procedure for a general barotropic 
equation of state and denoting by A the growth rate of 
the perturbation (such that Sp • 
eigenvalue equation 



e etc.), we obtain the 



d f p'jp) d 
dr \47rpr2 dr 



G_ 



A2 



Tlxir), 



(224) 



Anpr^ 

with appropriate boundary conditions (see Appendix EJ- 
This equation is equivalent to the Eddington (1926) equa- 
tion of pulsation but it proved to be more convenient in 
our previous investigations. We observe that Eqs. (|222|l 
and H224() are similar and that they coincide for marginal 
stability. This first dcmonstates the equivalence between 
the stability criterion of Eddington based on the equa- 
tion of pulsation ((224|l and the stability criterion of 
Chandrasekhar based on the minimization of the energy 
(|223|l . This also shows the equivalence between the maxi- 
mization of the J-function for a stellar system at fixed (3 
and AI and the dynamical stability (with respect to the 
Jeans-Euler equations) of the corresponding barotropic 
star. Note that we have only studied the dynamical sta- 
bility of barotropic stars with respect to small perturba- 
tions. However, we conjecture that the maximization of 
J at fixed M and (3 is a condition of nonlinear dynam- 
ical stability for the Jeans-Euler equation. The proof is 
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expected to be similar to the corresponding one for the 
Vlasov equation. 

7.9. A new/ interpretation of Antonov's first law 

We can now combine the previous results in order to pro- 
vide a new interpretation of Antonov's first law (Binney & 
Tremaine 1987) based on a thermodynamical analogy. We 
have seen in Sec. 17.71 that a stellar system with / = /(e) 
and /'(e) < is dynamically stable (nonlinearly) with re- 
spect to collisionless perturbations if / is a maximum of a 
7?- function S[f] at fixed E and M. According to the dis- 
cussion following Eqs. (|204() and 12051) . a sufficient con- 
dition of dynamical stability is that / is a maximum of 
the J-function J[f] at fixed (3 and M. Now, according 
to Sec. 17.81 this is equivalent to the dynamical stability 
of the corresponding barotropic gas. Therefore, we come 
to the conclusion that "a stellar system is dynamically 
stable whenever the corresponding barotropic star is sta- 
ble", which is Antonov's first law. Owing to the thermo- 
dynamical analogy, Antonov's first law is similar to the 
property that "canonical stability implies microcanonical 
stability" in thermodynamics. We note that the recipro- 
cal of Antonov's first law is not true in general. If a stellar 
system is stable this does not necessarily imply that the 
corresponding barotropic gas is stable. Similarly, in ther- 
modynamics, microcanonical stability does not necessarily 
imply canonical stability. This is the case only if the sta- 
tistical ensembles are equivalent. We know, however, that 
this equivalence is broken for self-gravitating systems. 

Let us illustrate these results in the case of polytropes. 
Taruya & Sakagami (2002a, b) have investigated the sta- 
bility of stellar systems in the framework of Tsallis gener- 
alized thermodynamics (see also Sec. l5.2f) . Due to the ther- 
modynamical analogy discussed in Sees. 17.71 and 17.81 the 
study of Taruya & Sakagami can be used to determine the 
dynamical stability of collisionless stellar polytropes and 
polytropic gaseous spheres. In the generalized thermody- 
namical approach, it is found that the statistical ensembles 
are inequivalent. This implies that the dynamical stabil- 
ity of stellar polytropes and gaseous polytropes will also 
be inequivalent (Chavanis 2002b). Self-gravitating systems 
described by Tsallis entropy are thermodynamically stable 
in the microcanonical ensemble if n < 5 or if the density 
contrast is sufficiently low for n > 5. In that case, they 
are maxima of Sq[f] at fixed M and E. According to Sec. 
17.71 this implies that stellar polytropes are dynamically 
stable solutions of the Vlasov equation for n < 5 (and if 
the density contrast is sufficiently low for n > 5). On the 
other hand, self-gravitating systems described by Tsallis 
entropy are thermodynamically stable in the canonical en- 
semble if n < 3 or if the density contrast is sufficiently low 
for n > 3. In that case, they are maxima of Jq[f] at fixed 
M and f3. According to Sec. 17.81 this implies that poly- 
tropic gaseous spheres are dynamically stable solutions 
of the Jeans-Euler equations if, and only, if n < 3 (or 
if, and only, if the density contrast is sufficiently low for 



n > 3). We note, in particular, that gaseous polytropes 
with 3 < ?i < 5 are unstable for sufficiently high density 
contrasts (e.g., complete polytropes) while corresponding 
stellar polytropes are stable. This shows that the recipro- 
cal of Antonov's first law is wrong for polytropes. Due to 
the thermodynamical analogy, this can be related to an 
inequivalence of statistical ensembles in a region of neg- 
ative specific heats (Chavanis 2002b). Our present study 
shows that the results obtained for isothermal stellar sys- 
tems and stellar polytropes can be generalized to any H- 
function. 

8. Conclusion 

We have performed an exhaustive study of the thermo- 
dynamics of self-gravitating systems in various ensembles. 
This paper completes previous investigations on the sub- 
ject and all ensembles have now been treated. Contrary to 
ordinary (extensive) systems, we have to perform a spe- 
cific study in each ensemble since the stability limits differ 
from one to the other. Remarkably, the thermodynami- 
cal stability problem can be studied analytically or with 
simple graphical constructions. The dynamical stability 
of isothermal gaseous spheres can also be studied analyt- 
ically both in Newtonian mechanics for the Euler- Jeans 
equations (Chavanis 2002a) and in general relativity for 
the Einstein equations (Chavanis 2002d). 

These results can be of relevance both for astrophysi- 
cists and statistical mechanicians and could make a bridge 
between these two communities. On an astrophysical point 
of view, they show that we must be careful to precisely 
define the ensemble in which we work since they are not 
equivalent. This does not affect the structure of the equi- 
librium configuration but it may affect its stability. On a 
physical point of view, this study fills an important gap 
in the statisical mechanics literature since the case of self- 
gravitating systems is not discussed at all in standard text- 
books of statistical mechanics and thermodynamics. 

We have also discussed the relevance of Tsallis general- 
ized thermodynamics for stellar systems. Collisionless stel- 
lar systems such as elliptical galaxies can achieve a metae- 
quilibrium state as a result of a violent relaxation. The 
Boltzmann entropy Ssif] is the correct entropy for these 
systems (in a coarse-grained sense and assuming that the 
system is non-degenerate) . It measures the disorder where 
the disorder is equal to the number of microstates consis- 
tent with a given macrostate. However, Ssif] is not maxi- 
mized by the system because of incomplete relaxation (this 
is independant on the fact that Ssif] has no maximum!). 
In any case, the state resulting from incomplete violent 
relaxation is a nonlinearly stable solution of the Vlasov 
equation (on a coarse-grained scale). If / = /(e), it max- 
imizes a H-function S[f] = — / C{f)d?'Y<fiw, where C(/) 
is a convex function, at fixed mass and energy. Therefore, 
we can use a thermodynamical analogy to study the dy- 
namical stability of collisionless stellar systems. Tsallis en- 
tropies Sq [/] are a particular class of H-functions leading 
to stellar polytropes. Stellar polytropes do not give a good 
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fit of elliptical galaxies. A better model of incomplete vio- 
lent relaxation, motivated by precise physical arguments, 
consists of an isothermal core and a polytropic halo with 
index n ~ 4 (Hjorth & Madsen 1993). This can be fitted 
by a composite Tsallis (polytropic) model with q ~ 1 in 
the core (complete mixing) and q ~ 7/5 in the envelope 
(incomplete mixing). More generally, the state resulting 
from incomplete violent relaxation does not necessarily 
maximize a H-function at fixed mass and energy. Due to 
the Jeans theorem, the distribution function can depend 
on other integrals than the stellar energy or the Jacobi 
energy. It is possible, nevertheless, that relevant distribu- 
tion functions arising from incomplete violent relaxation 
maximize a iJ-function at fixed mass, energy, angular mo- 
mentum and additional (adiabatic) invariants. 
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Appendix B: Some useful identities for polytropic 
spheres 

Using integrations by parts similar to those performed in 
Appendix we can easily establish the identities 
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Appendix A: Some useful identities for isothermal 
spheres 

Using the Eniden equation 114(1 . we have 

^^e-'^di = / {^'^ip'ydi ^ aV(a) = avo, (A.l) 

JO 

which establishes Eq. Ht)7|) . On the other hand, using an 
integration by parts, we obtain 

[%'ee~^di = - r e^(e-^)di 
Jo Jo d^ 



3e^e-'^de 



(A.2) 



Combining with Eq. (|A.1|I and introducing the Milne vari- 
ables, we establish Eq. (|68|l . To establish Eq. (|69() . we start 
from the identity 
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which results from a simple integration by parts. Thus, 
^2^'2d^ = -a^'iaf + 2 r ^'fe~^d^, (A.4) 



^0 

where we have used the Emden equation (|14|l . Combining 
this relation with Eq. ((A.2p , we obtain 



^^iP'^d^^avoi6-vo-2uo). 
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Now, 



e^e-^d^= / i^e^yd^ 
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From Eqs. I(A.5|I and l(A.6|l results Eq. 



Appendix C: Equation of pulsations and boundary 
conditions 

In this paper and previous ones, we have studied the dy- 
namical stability of barotropic gaseous spheres with the 
equation of pulsation 



dr \47r/9r-^ dr 
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Airpr- 



(C.l) 



where A is the growth rate of the perturbation (such 
that Sp 



Sv 



e-^* etc..) and q{r) = SM{r) 



SpAirr^dr represent the mass perturbation within the 
sphere of radius r which is related to the velocity by 
(Chavanis 2002a) 
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Obviously q{0) = and more precisely q ^ for r — s- 0. 
If the system is confined within a sphere of fixed radius 
R, the conservation of mass imposes q{R) = 0. This can 
also be obtained from Eq. l|C.2p if we set 6v{R) — 0. If the 
system is free and the density vanishes at its surface, i.e. 
p{R) — 0, like for a complete polytrope, then Eq. (|C.2|I 
again implies that q{R) — (Chavanis 2002b). If, finally, 
the gaseous sphere supports a fixed external pressure (iso- 
baric ensemble) the boundary condition requires that the 
Lagrangian derivative of the pressure vanishes at the sur- 
face of the configuration, i.e. 



d r r dp 
-jjOp + dv— = 0, 
dt dr 



at ?■ 



R. 



(C.3) 



Using dSp/dt = XSp and Eqs. HC.2|) and (|122|l . we get for 
an arbitrary barotropic equation of state 



dq _ qdp ^ ^ 
dr p dr ' 



at r = i?. 



(C.4) 



If we set Sv = ddr/dt = XSr and introduce the 
Lagrangian displacement ^(r) = 6r/r = 6v/Xr, we can 
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easily check that Eq. HC.ljl is equivalent to the Eddington 
(1926) equation of pulsations 



— \jrfr — 
ar \ ar 



ar 



4)p] 



(C.5) 



where j(r) = dlnp/dlnp. The boundary condition at the 
origin is ^'(0) = 0. On the other hand, if the system is 
enclosed within a box £,{R) — 0. Finally, if the pressure at 
the surface of the sphere is constant (possibly zero), the 
boundary condition 1C.4() becomes 



dr 



— ^ = — — , at r = i?. 



(C.6) 
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